theory-of-probability-and-mathematical-statistics/10. Критерии согласия, нормальности, однородности, независимости, случайности.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Гипотезы бывают: о виде распределения (критерии согласия), о однородности, о независимости, о случайности\n",
    "https://docs.scipy.org/doc/scipy/reference/stats.html#statistical-tests"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<strong>Автор:</strong> Маркашева Вера Алексеевна\\\n",
    "Косметические изменения: Царюк Артем Владимирович"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Критерии согласия: Колмогорова, Хи-квадрат, критерий квантилей и так далее\n",
    "\n",
    "Критерии нормальности: Жарки-Бера, Шапиро-Уилка, критерий квантилей реализуется на QQ-plot\n",
    "\n",
    "Критерии однородности: Смирнова, хи-квадрат\n",
    "\n",
    "Критерии независимости: Хи-квадрат, Спирмена, Кендалла\n",
    "\n",
    "Критерии случайности: критерий инверсий"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#  Простые гипотезы "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Умение тестировать статистические гипотезы – один из наиболее востребованных навыков аналитика. Давайте приведем пример. Мы сделали новую модель, которая умеет предсказывать, совершит ли покупку на сайте человек в ближайшую неделю. Как проверить, что мы хорошо научились угадывать? Очень просто - давайте для наших пользователей сгенерируем предсказания модели, а потом сравним эти предсказания с фактом покупок. Логично. Наша модель будет полезна, если она работает лучше, чем случайная угадайка (1/2, что купит и 1/2, что не купит). \n",
    "\n",
    "Давайте проверим гипотезу, что наша модель == подбрасывание монетки. Говорят данные в пользу такого утверждения или против него?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Проверка гипотезы формализуется так:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Выборка**   $X_n = {X_1 ... X_n}$ из неизвестного распределения $X \\sim P$\n",
    "\n",
    "**Нулевая гипотеза**    $H_0: P \\in \\omega$ ,  где $\\omega$ - известное семейство распределений\n",
    "\n",
    "**Альтернативная гипотеза**     $H_1: P \\notin \\omega$ \n",
    "\n",
    "**Статистика**    $T(X_n)$\n",
    "\n",
    "**Нулевое распределение**    $F(T)$ - распределение статистики в случае справедливости $H_0$\n",
    "\n",
    "Если справедлива гипотеза $H_1$ , то распределение нашей статистики имеет ДРУГОЕ распределение, не  $F(T)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Определение** *p-value* или достигаемый уровень значимости это вероятность получить такое же или еще более экстремальное значение статистики $T(X_n)$, чем было получено в эксперименте, при условии справедливости $H_0$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "То есть если p-value достаточно мало, это значит, что данные свидетельствуют ПРОТИВ нулевой гипотезы, потому что вероятность получить такое значение при ее справедливости – маленькая величина. Обычно p-value сравнивают с наперед заданным уровнем значимости $\\alpha$, которы обычно берется за 0.05 или 0.01. Ну и если $p-value < \\alpha$, говорят, что нулевая гипотеза ОТВЕРГАЕТСЯ в пользу альтернативы."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "При тестировании гипотез существуют два вида ошибок: неверно отклоняется или неверно принимается гипотеза. В тестировании гипотез считается, что ошибка первого рода критичнее, так как вероятность отвержения нулевой гипотезы оказывается ограниченной сверху нашим уровнем значимости:\n",
    "$$P(\\text{ошибка первого рода})\\le\\alpha$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ошибки второго рода так явно не минимизируются. Вероятность ошибки второго рода напрямую связана с понятием мощности статистического критерия,  то есть:\n",
    "$$P(\\text{вероятность отвергнуть $H_0$}|H_1) = 1 - P(\\text{вероятность принять $H_0$}|H_1)$$\n",
    "\n",
    "То есть среди всех подходящих критериев мы должны выбрать тот, у которого максимальная мощность."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Даже в терминологии здесь видна несимметричность. Если достигаемый уровень значимости меньше $\\alpha$, то гипотеза $H_0$ отвергается, а если не меньше, то НЕ отвергается (а не принимается, как ожидалось). То есть это не доказательство того, что $H_0$ верна."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Размер эффекта"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Когда мы проводим тестирование гипотез, нас на самом деле интересует некий дополнительный эффект, который мы можем получить за счет тестируемой функциональности, например кол-во дополнительных заказов от новой программы лояльности, увеличение длины диалога с голосовым помощником от новой диалоговой механики или добавочная вероятность выздоровления от применения нового способа лечения болезни. Мы с вами не можем оценить этот эффект на всей генеральной совокупности и поэтому оцениваем его по выборке. Вот тут-то нам и пригождается достигаемый уровень значимости, чтобы понять, с какой вероятностью такой же эффект можно было бы получить случайно."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "p сильно зависит не только от размера эффекта, но и от размера выборки. На больших выборках вы будете замечать все более и более тонкие отклонения от нулевой гипотезы и при каком-то n почти все $H_0$ будут отвергаться."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Поэтому эффект часто надо оценивать с прикладной точки зрения. Например если в течение трех лет регулярные тренировки в спортзале по некой системе помогли людям с ожирением снизить массу тела в среднем на 200 граммов и этот эффект статистически значим, представляется весьма сомнительным, что кто-либо желающий похудеть получит от такой цифры моральное удовлетворение. А значит в контексте похудения этот эффект практически бесполезен."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "С другой стороны, не все статистически не значимые эксперименты бессмысленны. Скажем, если вы изобрели новое лекарство от вируса COVID-19, протестировали его на добровольцах, оно показывает хорошие результаты (люди выздоравливают, побочных эффектов нет), но при этом результат не стал статистически значимым, это не повод все бросить. Исследования стоит продолжать, потому что практически такой эффект может оказаться значимым. Возможно вам стоит подумать об уменьшении дисперсии в вашем исследовании."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "В любом случае, статистически значимый результат дает нечто, на что можно опираться, как на наиболее вероятное событие в будущем.\n",
    "Уровень значимости стоит понимать как минимум как меру доверия к возможному событию."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Гипотеза о виде распределения\n",
    "\n",
    "# Параметрические критерии \n",
    "Так называют критерии, которые делают некие предположения о характере распределений. А мы, в частности, будем оценивать какие-либо параметры этих распределений. Начнем с проверки гипотез, связанных с признаками, распределенными нормально. Начнем с так называемых критериев Стьюдента о математических ожиданиях. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Биномиальный тест для доли"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Lady tasting tea experiment\n",
    "\n",
    "Известный статистик Рональд Фишер провел эксперимент, в котором предложил даме 8 чашек чая, в 4 из которых сначала добавили чай, а потом молоко, а в 4 другие – наоборот (сначала молоко, потом чай). Ей необходимо было разбить чашки на две группы (первым было молоко или первым был чай). Проверялась гипотеза о том, что дама случайно угадывает, что было налито первым. В ходе анализа этого эксперимента Фишер придумал так называемый [Точный тест Фишера](https://ru.wikipedia.org/wiki/Точный_тест_Фишера), а дама подтвердила свое умение определять, что было первым налито в чашку.\n",
    "\n",
    "Проведем эксперимент чуть иначе. Предложим даме 20 чашек (по две за раз) и будем спрашивать в какой чашке из пары сначала налито молоко. Тогда наш эксперимент - бинарный вектор длины 10 (допустим 0 - налит чай, а 1 -налито молоко)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Если нулевая гипотеза о выборе наугад справедлива , то мы равновероятно можем получить любой из $2^n$ бинарных векторов длины $n$. \n",
    "\n",
    "Нулевая гипотеза выбора наугад соответствует значению $p=\\frac1{2}$, то есть, в каждом эксперименте вероятность выбора чая, в который сначала налито молоко – $\\frac1{2}$. Сумма $n$ одинаково распределённых бернуллиевских случайных величин с параметром $p$ имеет биномиальное распределение $Bin(n, p)$. Значит наше нулевое распределение статистики $T$ — $Bin\\left(n, \\frac1{2}\\right)$.\n",
    "\n",
    "Мы уже сказали, что $n=10.$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from matplotlib import pylab\n",
    "from scipy import stats\n",
    "import scipy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 10\n",
    "F_H0 = stats.binom(n, 0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiwAAAGdCAYAAAAxCSikAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAiZElEQVR4nO3df1DVVeL/8Rc/4uJPUkkuGAamZZqAiTKYfZwd73hxnJKtDB0Lo8ZmXdm0u0tGK+Bmu6CZSyYjq7vu6kymNbO624+l7K6424SyQa5rWWmrodK9qLuK4ggO3O8fjdfvXfHHRfAeL8/HzJ2VN+eee86gt+e+ed97Qzwej0cAAAAGCw30AgAAAK6GYAEAAMYjWAAAgPEIFgAAYDyCBQAAGI9gAQAAxiNYAACA8QgWAABgvPBAL6AztLW1qb6+Xn369FFISEiglwMAAK6Bx+PR6dOnFRcXp9DQK59DCYpgqa+vV3x8fKCXAQAAOuDw4cO6/fbbrzgmKIKlT58+kr7fcN++fQO8GgAAcC0aGxsVHx/v/e/4lQRFsFz4NVDfvn0JFgAAbjLXcjkHF90CAADjESwAAMB4BAsAADAewQIAAIxHsAAAAON1KFjKysqUkJCgyMhIpaWlqbq6+rJj165dqwceeED9+vVTv379ZLPZLhn/5JNPKiQkxOeWkZHRkaUBAIAg5HewbN68WQ6HQ0VFRaqtrVVycrLsdrsaGhraHV9ZWamZM2dq+/btqqqqUnx8vCZPnqyjR4/6jMvIyNB3333nvb355psd2xEAAAg6IR6Px+PPHdLS0jR27FitWrVK0vdvix8fH6+f/OQneuGFF656/9bWVvXr10+rVq1Sdna2pO/PsJw8eVJbt271fwf6/o1noqKidOrUKd6HBQCAm4Q///326wxLS0uLampqZLPZLk4QGiqbzaaqqqprmuPs2bM6f/68+vfv73O8srJSAwcO1N133625c+fqxIkTl52jublZjY2NPjcAABC8/AqW48ePq7W1VTExMT7HY2Ji5HK5rmmOhQsXKi4uzid6MjIytGHDBjmdTi1dulQ7duzQlClT1Nra2u4cxcXFioqK8t74HCEAAILbDX1r/pKSEm3atEmVlZWKjIz0Hp8xY4b3z6NGjVJSUpLuvPNOVVZWatKkSZfMk5+fL4fD4f36wmcRAACA4OTXGZbo6GiFhYXJ7Xb7HHe73bJarVe87/Lly1VSUqIPP/xQSUlJVxw7ZMgQRUdH68CBA+1+32KxeD83iM8PAgAg+PkVLBERERozZoycTqf3WFtbm5xOp9LT0y97v2XLlmnJkiWqqKhQamrqVR/nyJEjOnHihGJjY/1ZHgAACFJ+v6zZ4XBo7dq1Wr9+vfbt26e5c+eqqalJOTk5kqTs7Gzl5+d7xy9dulQFBQVat26dEhIS5HK55HK5dObMGUnSmTNnlJeXp507d+rQoUNyOp2aNm2ahg4dKrvd3knbBAAANzO/r2HJysrSsWPHVFhYKJfLpZSUFFVUVHgvxK2rq1No6MUOWr16tVpaWvToo4/6zFNUVKTFixcrLCxMe/bs0fr163Xy5EnFxcVp8uTJWrJkiSwWy3VuD+i4hBfe69L5D5VM7dL5ASCY+P0+LCbifVjQFQgWAOhaXfY+LAAAAIFAsAAAAOMRLAAAwHgECwAAMB7BAgAAjEewAAAA4xEsAADAeAQLAAAwHsECAACMR7AAAADjESwAAMB4BAsAADAewQIAAIxHsAAAAOMRLAAAwHgECwAAMB7BAgAAjEewAAAA4xEsAADAeAQLAAAwHsECAACMR7AAAADjESwAAMB4BAsAADAewQIAAIxHsAAAAOMRLAAAwHgECwAAMB7BAgAAjEewAAAA4xEsAADAeAQLAAAwHsECAACMR7AAAADjESwAAMB4BAsAADAewQIAAIxHsAAAAOMRLAAAwHgECwAAMB7BAgAAjEewAAAA4xEsAADAeAQLAAAwHsECAACMR7AAAADjESwAAMB4BAsAADAewQIAAIxHsAAAAOMRLAAAwHgECwAAMB7BAgAAjEewAAAA4xEsAADAeAQLAAAwHsECAACMFx7oBQDXIuGF97ps7kMlU7tsbn915T4ls/YKAP7gDAsAADAewQIAAIzXoWApKytTQkKCIiMjlZaWpurq6suOXbt2rR544AH169dP/fr1k81mu2S8x+NRYWGhYmNj1aNHD9lsNu3fv78jSwMAAEHI72DZvHmzHA6HioqKVFtbq+TkZNntdjU0NLQ7vrKyUjNnztT27dtVVVWl+Ph4TZ48WUePHvWOWbZsmVauXKny8nLt2rVLvXr1kt1u17lz5zq+MwAAEDT8DpYVK1Zozpw5ysnJ0YgRI1ReXq6ePXtq3bp17Y5/44039OMf/1gpKSkaPny4fvvb36qtrU1Op1PS92dXSktLtWjRIk2bNk1JSUnasGGD6uvrtXXr1uvaHAAACA5+BUtLS4tqampks9kuThAaKpvNpqqqqmua4+zZszp//rz69+8vSTp48KBcLpfPnFFRUUpLS7vsnM3NzWpsbPS5AQCA4OVXsBw/flytra2KiYnxOR4TEyOXy3VNcyxcuFBxcXHeQLlwP3/mLC4uVlRUlPcWHx/vzzYAAMBN5oa+SqikpESbNm3Sli1bFBkZ2eF58vPzderUKe/t8OHDnbhKAABgGr/eOC46OlphYWFyu90+x91ut6xW6xXvu3z5cpWUlOijjz5SUlKS9/iF+7ndbsXGxvrMmZKS0u5cFotFFovFn6UDAICbmF9nWCIiIjRmzBjvBbOSvBfQpqenX/Z+y5Yt05IlS1RRUaHU1FSf7yUmJspqtfrM2djYqF27dl1xTgAA0H34/db8DodDs2fPVmpqqsaNG6fS0lI1NTUpJydHkpSdna1BgwapuLhYkrR06VIVFhZq48aNSkhI8F6X0rt3b/Xu3VshISFasGCBXn75ZQ0bNkyJiYkqKChQXFycMjMzO2+nAADgpuV3sGRlZenYsWMqLCyUy+VSSkqKKioqvBfN1tXVKTT04omb1atXq6WlRY8++qjPPEVFRVq8eLEk6fnnn1dTU5OeeeYZnTx5UhMmTFBFRcV1XecCAACCR4c+/DA3N1e5ubntfq+ystLn60OHDl11vpCQEL300kt66aWXOrIcAAAQ5PgsIQAAYDyCBQAAGI9gAQAAxiNYAACA8QgWAABgPIIFAAAYj2ABAADGI1gAAIDxCBYAAGA8ggUAABiPYAEAAMYjWAAAgPEIFgAAYDyCBQAAGI9gAQAAxiNYAACA8QgWAABgPIIFAAAYj2ABAADGI1gAAIDxCBYAAGA8ggUAABiPYAEAAMYjWAAAgPEIFgAAYDyCBQAAGI9gAQAAxiNYAACA8QgWAABgPIIFAAAYj2ABAADGI1gAAIDxCBYAAGA8ggUAABiPYAEAAMYjWAAAgPEIFgAAYDyCBQAAGI9gAQAAxiNYAACA8QgWAABgPIIFAAAYj2ABAADGI1gAAIDxCBYAAGA8ggUAABiPYAEAAMYjWAAAgPEIFgAAYDyCBQAAGI9gAQAAxiNYAACA8QgWAABgPIIFAAAYj2ABAADGI1gAAIDxCBYAAGA8ggUAABiPYAEAAMYjWAAAgPEIFgAAYLwOBUtZWZkSEhIUGRmptLQ0VVdXX3bs559/rkceeUQJCQkKCQlRaWnpJWMWL16skJAQn9vw4cM7sjQAABCE/A6WzZs3y+FwqKioSLW1tUpOTpbdbldDQ0O748+ePashQ4aopKREVqv1svOOHDlS3333nff28ccf+7s0AAAQpPwOlhUrVmjOnDnKycnRiBEjVF5erp49e2rdunXtjh87dqxeeeUVzZgxQxaL5bLzhoeHy2q1em/R0dH+Lg0AAAQpv4KlpaVFNTU1stlsFycIDZXNZlNVVdV1LWT//v2Ki4vTkCFDNGvWLNXV1V12bHNzsxobG31uAAAgePkVLMePH1dra6tiYmJ8jsfExMjlcnV4EWlpafrDH/6giooKrV69WgcPHtQDDzyg06dPtzu+uLhYUVFR3lt8fHyHHxsAAJjPiFcJTZkyRdOnT1dSUpLsdrvef/99nTx5Um+99Va74/Pz83Xq1Cnv7fDhwzd4xQAA4EYK92dwdHS0wsLC5Ha7fY673e4rXlDrr1tvvVV33XWXDhw40O73LRbLFa+HAQAAwcWvMywREREaM2aMnE6n91hbW5ucTqfS09M7bVFnzpzRN998o9jY2E6bEwAA3Lz8OsMiSQ6HQ7Nnz1ZqaqrGjRun0tJSNTU1KScnR5KUnZ2tQYMGqbi4WNL3F+p+8cUX3j8fPXpUu3fvVu/evTV06FBJ0s9+9jM9+OCDuuOOO1RfX6+ioiKFhYVp5syZnbVPAABwE/M7WLKysnTs2DEVFhbK5XIpJSVFFRUV3gtx6+rqFBp68cRNfX29Ro8e7f16+fLlWr58uSZOnKjKykpJ0pEjRzRz5kydOHFCt912myZMmKCdO3fqtttuu87tAQCAYOB3sEhSbm6ucnNz2/3ehQi5ICEhQR6P54rzbdq0qSPLAAAA3YQRrxICAAC4EoIFAAAYj2ABAADGI1gAAIDxCBYAAGA8ggUAABiPYAEAAMYjWAAAgPEIFgAAYDyCBQAAGI9gAQAAxiNYAACA8QgWAABgPIIFAAAYj2ABAADGI1gAAIDxCBYAAGA8ggUAABiPYAEAAMYjWAAAgPEIFgAAYDyCBQAAGI9gAQAAxiNYAACA8QgWAABgPIIFAAAYj2ABAADGI1gAAIDxCBYAAGA8ggUAABiPYAEAAMYjWAAAgPEIFgAAYDyCBQAAGI9gAQAAxiNYAACA8QgWAABgPIIFAAAYj2ABAADGI1gAAIDxwgO9ANx8El54r8vmPlQytcvmxuXxMwVgOs6wAAAA4xEsAADAeAQLAAAwHsECAACMR7AAAADjESwAAMB4BAsAADAewQIAAIxHsAAAAOMRLAAAwHgECwAAMB7BAgAAjEewAAAA4xEsAADAeAQLAAAwHsECAACMR7AAAADjESwAAMB4BAsAADBeh4KlrKxMCQkJioyMVFpamqqrqy879vPPP9cjjzyihIQEhYSEqLS09LrnBAAA3YvfwbJ582Y5HA4VFRWptrZWycnJstvtamhoaHf82bNnNWTIEJWUlMhqtXbKnAAAoHvxO1hWrFihOXPmKCcnRyNGjFB5ebl69uypdevWtTt+7NixeuWVVzRjxgxZLJZOmRMAAHQvfgVLS0uLampqZLPZLk4QGiqbzaaqqqoOLaAjczY3N6uxsdHnBgAAgpdfwXL8+HG1trYqJibG53hMTIxcLleHFtCROYuLixUVFeW9xcfHd+ixAQDAzeGmfJVQfn6+Tp065b0dPnw40EsCAABdKNyfwdHR0QoLC5Pb7fY57na7L3tBbVfMabFYLns9DAAACD5+nWGJiIjQmDFj5HQ6vcfa2trkdDqVnp7eoQV0xZwAACC4+HWGRZIcDodmz56t1NRUjRs3TqWlpWpqalJOTo4kKTs7W4MGDVJxcbGk7y+q/eKLL7x/Pnr0qHbv3q3evXtr6NCh1zQnAADo3vwOlqysLB07dkyFhYVyuVxKSUlRRUWF96LZuro6hYZePHFTX1+v0aNHe79evny5li9frokTJ6qysvKa5gQAAN2b38EiSbm5ucrNzW33exci5IKEhAR5PJ7rmhMAAHRvN+WrhAAAQPdCsAAAAOMRLAAAwHgECwAAMB7BAgAAjEewAAAA4xEsAADAeAQLAAAwHsECAACMR7AAAADjESwAAMB4BAsAADAewQIAAIxHsAAAAOMRLAAAwHgECwAAMB7BAgAAjEewAAAA4xEsAADAeAQLAAAwHsECAACMR7AAAADjESwAAMB4BAsAADAewQIAAIxHsAAAAOMRLAAAwHgECwAAMB7BAgAAjEewAAAA4xEsAADAeAQLAAAwHsECAACMR7AAAADjESwAAMB4BAsAADAewQIAAIxHsAAAAOMRLAAAwHgECwAAMB7BAgAAjEewAAAA4xEsAADAeAQLAAAwHsECAACMR7AAAADjESwAAMB4BAsAADAewQIAAIxHsAAAAOMRLAAAwHgECwAAMB7BAgAAjEewAAAA4xEsAADAeAQLAAAwHsECAACMR7AAAADjESwAAMB4BAsAADAewQIAAIzXoWApKytTQkKCIiMjlZaWpurq6iuOf/vttzV8+HBFRkZq1KhRev/9932+/+STTyokJMTnlpGR0ZGlAQCAIOR3sGzevFkOh0NFRUWqra1VcnKy7Ha7Ghoa2h3/ySefaObMmXr66af12WefKTMzU5mZmdq7d6/PuIyMDH333Xfe25tvvtmxHQEAgKDjd7CsWLFCc+bMUU5OjkaMGKHy8nL17NlT69ata3f8a6+9poyMDOXl5emee+7RkiVLdN9992nVqlU+4ywWi6xWq/fWr1+/ju0IAAAEHb+CpaWlRTU1NbLZbBcnCA2VzWZTVVVVu/epqqryGS9Jdrv9kvGVlZUaOHCg7r77bs2dO1cnTpy47Dqam5vV2NjocwMAAMHLr2A5fvy4WltbFRMT43M8JiZGLper3fu4XK6rjs/IyNCGDRvkdDq1dOlS7dixQ1OmTFFra2u7cxYXFysqKsp7i4+P92cbAADgJhMe6AVI0owZM7x/HjVqlJKSknTnnXeqsrJSkyZNumR8fn6+HA6H9+vGxkaiBQCAIObXGZbo6GiFhYXJ7Xb7HHe73bJare3ex2q1+jVekoYMGaLo6GgdOHCg3e9bLBb17dvX5wYAAIKXX8ESERGhMWPGyOl0eo+1tbXJ6XQqPT293fukp6f7jJekbdu2XXa8JB05ckQnTpxQbGysP8sDAABByu9XCTkcDq1du1br16/Xvn37NHfuXDU1NSknJ0eSlJ2drfz8fO/4+fPnq6KiQq+++qq+/PJLLV68WJ9++qlyc3MlSWfOnFFeXp527typQ4cOyel0atq0aRo6dKjsdnsnbRMAANzM/L6GJSsrS8eOHVNhYaFcLpdSUlJUUVHhvbC2rq5OoaEXO2j8+PHauHGjFi1apBdffFHDhg3T1q1bde+990qSwsLCtGfPHq1fv14nT55UXFycJk+erCVLlshisXTSNgEAwM2sQxfd5ubmes+Q/K/KyspLjk2fPl3Tp09vd3yPHj30wQcfdGQZAACgm+CzhAAAgPGMeFkzOi7hhfe6bO5DJVO7bG6Av7sA/MEZFgAAYDyCBQAAGI9gAQAAxiNYAACA8QgWAABgPIIFAAAYj2ABAADGI1gAAIDxCBYAAGA8ggUAABiPYAEAAMYjWAAAgPEIFgAAYDyCBQAAGI9gAQAAxiNYAACA8QgWAABgPIIFAAAYj2ABAADGI1gAAIDxCBYAAGA8ggUAABiPYAEAAMYjWAAAgPEIFgAAYDyCBQAAGI9gAQAAxiNYAACA8QgWAABgPIIFAAAYj2ABAADGI1gAAIDxCBYAAGA8ggUAABiPYAEAAMYjWAAAgPEIFgAAYDyCBQAAGI9gAQAAxiNYAACA8QgWAABgPIIFAAAYj2ABAADGI1gAAIDxCBYAAGA8ggUAABiPYAEAAMYjWAAAgPHCA72AYJLwwntdNvehkqldNjfQXfBvFLh5cYYFAAAYj2ABAADGI1gAAIDxCBYAAGA8ggUAABiPYAEAAMYjWAAAgPEIFgAAYDyCBQAAGK9DwVJWVqaEhARFRkYqLS1N1dXVVxz/9ttva/jw4YqMjNSoUaP0/vvv+3zf4/GosLBQsbGx6tGjh2w2m/bv39+RpQEAgCDkd7Bs3rxZDodDRUVFqq2tVXJysux2uxoaGtod/8knn2jmzJl6+umn9dlnnykzM1OZmZnau3evd8yyZcu0cuVKlZeXa9euXerVq5fsdrvOnTvX8Z0BAICg4XewrFixQnPmzFFOTo5GjBih8vJy9ezZU+vWrWt3/GuvvaaMjAzl5eXpnnvu0ZIlS3Tfffdp1apVkr4/u1JaWqpFixZp2rRpSkpK0oYNG1RfX6+tW7de1+YAAEBw8OvDD1taWlRTU6P8/HzvsdDQUNlsNlVVVbV7n6qqKjkcDp9jdrvdGyMHDx6Uy+WSzWbzfj8qKkppaWmqqqrSjBkzLpmzublZzc3N3q9PnTolSWpsbJQk3Vv0gT/b8tveX9jbPd7WfLbLHvPC3njM4HjMyz1uIB6zqx+3uz8mgMu78O/G4/FcdaxfwXL8+HG1trYqJibG53hMTIy+/PLLdu/jcrnaHe9yubzfv3DscmP+V3FxsX7xi19ccjw+Pv7aNnKdokpvyMPwmEH8mIF6XB4zuB4TCBanT59WVFTUFcf4FSymyM/P9zlr09bWpv/85z8aMGCAQkJC/JqrsbFR8fHxOnz4sPr27dvZSzVGd9mn1H32yj6DT3fZa3fZp9R99trRfXo8Hp0+fVpxcXFXHetXsERHRyssLExut9vnuNvtltVqbfc+Vqv1iuMv/K/b7VZsbKzPmJSUlHbntFgsslgsPsduvfVWf7Zyib59+wb1X6YLuss+pe6zV/YZfLrLXrvLPqXus9eO7PNqZ1Yu8Oui24iICI0ZM0ZOp9N7rK2tTU6nU+np6e3eJz093We8JG3bts07PjExUVar1WdMY2Ojdu3addk5AQBA9+L3r4QcDodmz56t1NRUjRs3TqWlpWpqalJOTo4kKTs7W4MGDVJxcbEkaf78+Zo4caJeffVVTZ06VZs2bdKnn36qNWvWSJJCQkK0YMECvfzyyxo2bJgSExNVUFCguLg4ZWZmdt5OAQDATcvvYMnKytKxY8dUWFgol8ullJQUVVRUeC+araurU2joxRM348eP18aNG7Vo0SK9+OKLGjZsmLZu3ap7773XO+b5559XU1OTnnnmGZ08eVITJkxQRUWFIiMjO2GLV2axWFRUVHTJr5iCTXfZp9R99so+g0932Wt32afUffZ6I/YZ4rmW1xIBAAAEEJ8lBAAAjEewAAAA4xEsAADAeAQLAAAwXrcPlrKyMiUkJCgyMlJpaWmqrq4O9JI6VXFxscaOHas+ffpo4MCByszM1FdffRXoZXW5kpIS70vmg9HRo0f1+OOPa8CAAerRo4dGjRqlTz/9NNDL6lStra0qKChQYmKievTooTvvvFNLliy5ps8cMd3f/vY3Pfjgg4qLi1NISMglH/Tq8XhUWFio2NhY9ejRQzabTfv37w/MYq/DlfZ5/vx5LVy4UKNGjVKvXr0UFxen7Oxs1dfXB27BHXS1n+f/70c/+pFCQkJUWlp6w9bXma5lr/v27dNDDz2kqKgo9erVS2PHjlVdXd11P3a3DpbNmzfL4XCoqKhItbW1Sk5Olt1uV0NDQ6CX1ml27NihefPmaefOndq2bZvOnz+vyZMnq6mpKdBL6zL/+Mc/9Jvf/EZJSUmBXkqX+O9//6v7779ft9xyi/7yl7/oiy++0Kuvvqp+/foFemmdaunSpVq9erVWrVqlffv2aenSpVq2bJlef/31QC/tujU1NSk5OVllZWXtfn/ZsmVauXKlysvLtWvXLvXq1Ut2u13nzp27wSu9Plfa59mzZ1VbW6uCggLV1tbqj3/8o7766is99NBDAVjp9bnaz/OCLVu2aOfOndf0NvSmutpev/nmG02YMEHDhw9XZWWl9uzZo4KCgs55mxJPNzZu3DjPvHnzvF+3trZ64uLiPMXFxQFcVddqaGjwSPLs2LEj0EvpEqdPn/YMGzbMs23bNs/EiRM98+fPD/SSOt3ChQs9EyZMCPQyutzUqVM9Tz31lM+xhx9+2DNr1qwArahrSPJs2bLF+3VbW5vHarV6XnnlFe+xkydPeiwWi+fNN98MwAo7x//usz3V1dUeSZ5vv/32xiyqC1xun0eOHPEMGjTIs3fvXs8dd9zh+fWvf33D19bZ2ttrVlaW5/HHH++Sx+u2Z1haWlpUU1Mjm83mPRYaGiqbzaaqqqoArqxrnTp1SpLUv3//AK+ka8ybN09Tp071+bkGmz//+c9KTU3V9OnTNXDgQI0ePVpr164N9LI63fjx4+V0OvX1119Lkv75z3/q448/1pQpUwK8sq518OBBuVwun7/DUVFRSktLC+rnJun756eQkJDr/mw407S1temJJ55QXl6eRo4cGejldJm2tja99957uuuuu2S32zVw4EClpaVd8Vdk/ui2wXL8+HG1trZ636H3gpiYGLlcrgCtqmu1tbVpwYIFuv/++33eaThYbNq0SbW1td6PhQhW//73v7V69WoNGzZMH3zwgebOnatnn31W69evD/TSOtULL7ygGTNmaPjw4brllls0evRoLViwQLNmzQr00rrUheef7vTcJEnnzp3TwoULNXPmzKD7kMClS5cqPDxczz77bKCX0qUaGhp05swZlZSUKCMjQx9++KF++MMf6uGHH9aOHTuue36/35ofN6958+Zp7969+vjjjwO9lE53+PBhzZ8/X9u2bbshH+kQSG1tbUpNTdWvfvUrSdLo0aO1d+9elZeXa/bs2QFeXed566239MYbb2jjxo0aOXKkdu/erQULFiguLi6o9onvL8B97LHH5PF4tHr16kAvp1PV1NTotddeU21trUJCQgK9nC7V1tYmSZo2bZqee+45SVJKSoo++eQTlZeXa+LEidc1f7c9wxIdHa2wsDC53W6f4263W1arNUCr6jq5ubl69913tX37dt1+++2BXk6nq6mpUUNDg+677z6Fh4crPDxcO3bs0MqVKxUeHq7W1tZAL7HTxMbGasSIET7H7rnnnk65Ct8keXl53rMso0aN0hNPPKHnnnsu6M+gXXj+6S7PTRdi5dtvv9W2bduC7uzK3//+dzU0NGjw4MHe56Zvv/1WP/3pT5WQkBDo5XWq6OhohYeHd9nzU7cNloiICI0ZM0ZOp9N7rK2tTU6nU+np6QFcWefyeDzKzc3Vli1b9Ne//lWJiYmBXlKXmDRpkv71r39p9+7d3ltqaqpmzZql3bt3KywsLNBL7DT333//JS9N//rrr3XHHXcEaEVd4+zZsz4fpCpJYWFh3v8XF6wSExNltVp9npsaGxu1a9euoHpuki7Gyv79+/XRRx9pwIABgV5Sp3viiSe0Z88en+emuLg45eXl6YMPPgj08jpVRESExo4d22XPT936V0IOh0OzZ89Wamqqxo0bp9LSUjU1NSknJyfQS+s08+bN08aNG/WnP/1Jffr08f4OPCoqSj169Ajw6jpPnz59Lrkup1evXhowYEDQXa/z3HPPafz48frVr36lxx57TNXV1VqzZo3WrFkT6KV1qgcffFC//OUvNXjwYI0cOVKfffaZVqxYoaeeeirQS7tuZ86c0YEDB7xfHzx4ULt371b//v01ePBgLViwQC+//LKGDRumxMREFRQUKC4uTpmZmYFbdAdcaZ+xsbF69NFHVVtbq3fffVetra3e56f+/fsrIiIiUMv229V+nv8bYrfccousVqvuvvvuG73U63a1vebl5SkrK0v/93//px/84AeqqKjQO++8o8rKyut/8C557dFN5PXXX/cMHjzYExER4Rk3bpxn586dgV5Sp5LU7u33v/99oJfW5YL1Zc0ej8fzzjvveO69916PxWLxDB8+3LNmzZpAL6nTNTY2eubPn+8ZPHiwJzIy0jNkyBDPz3/+c09zc3Ogl3bdtm/f3u6/y9mzZ3s8nu9f2lxQUOCJiYnxWCwWz6RJkzxfffVVYBfdAVfa58GDBy/7/LR9+/ZAL90vV/t5/q+b+WXN17LX3/3ud56hQ4d6IiMjPcnJyZ6tW7d2ymOHeDxB8LaRAAAgqHXba1gAAMDNg2ABAADGI1gAAIDxCBYAAGA8ggUAABiPYAEAAMYjWAAAgPEIFgAAYDyCBQAAGI9gAQAAxiNYAACA8QgWAABgvP8HkGLZ0tHgEagAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.linspace(0,9,10)\n",
    "\n",
    "pylab.bar(x, F_H0.pmf(x))\n",
    "pylab.xlim(-0.5, 16.5) \n",
    "pylab.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Односторонняя альтернатива\n",
    "\n",
    "Означает, что мы проверяем гипотезу, что дама угадывает лучше чем случайный классификатор (чаще действительно угадывает молоко).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "BinomTestResult(k=9, n=10, alternative='greater', statistic=0.9, pvalue=0.0107421875)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stats.binomtest(9, 10, 0.5, alternative = 'greater')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Двусторонняя альтернатива\n",
    "\n",
    "Означает, что мы проверяем гипотезу, что угадывает не как случайный классификатор (или чаще действительно угадывает молоко или чаще вместо молока указывает на чай)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "BinomTestResult(k=9, n=10, alternative='two-sided', statistic=0.9, pvalue=0.021484375)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stats.binomtest(9, 10, 0.5, alternative = 'two-sided')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Заметьте, что и в первом случае для стат.значимости надо угадать 9 раз из 10, а во втором - или угадать 9 раз или ошибиться 9 раз."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Связь с доверительными интервалами"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Понятно, что тестирование гипотез напрямую связано с построением доверительных интервалов. Возьмем задачку о сравнении классификатора на сбалансированной выборке(50 на 50) с подбрасыванием монетки (генератором случайных чисел). Если построить 95% доверительный интервал для доли  верно угаданных ответов и он не будет содержать 0.5 (вероятность случайного угадывания), можно смело утверждать, что на уровне 95% гипотеза об идентичности вашего классификатора и генератора случайных чисел отвергается. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Выборка**    $X_n = {X_1 ... X_n}$ из неизвестного распределения $X \\sim P$\n",
    "\n",
    "**Нулевая гипотеза**    $H_0: \\theta = \\theta_0$\n",
    "\n",
    "**Альтернативная гипотеза**     $H_1: \\theta \\not= \\theta_0$ \n",
    "\n",
    "\n",
    "$H_0$ отвергается на уровне значимости $100(1-\\alpha)$, если доверительный интервал для $\\theta$ не содержит $\\theta_0$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*p-value* здесь можно найти численно как такое значение $\\alpha$, что $\\theta_0$ будет попадать ровно на границу доверительного интервала."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Пусть у нас есть ответы классификатора для задачи с двумя сбалансированными выборками. Будем обозначать правильный ответ за 1, а неправильный - за 0. На тестовой выборке длины 1000 классификатор верно угадал 600 из 1000 ответов. Лучше ли он, чем бросание монетки?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "BinomTestResult(k=600, n=1000, alternative='greater', statistic=0.6, pvalue=1.364232078032988e-10)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stats.binomtest(600, 1000, 0.5, alternative = 'greater')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Как проверить нормальность (и вообще принадлежность к распределению)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Мы до сих пор работали в предположении, что наша выборка взята из нормального распределения. На самом деле, это утверждение неплохо было бы проверять. Например это можно сделать с помощью критерия Хи-квадрат. Здесь мы его применим для сравнения с нормальным распределением, но вообще он годится для проверки соответствия распределения ЛЮБОМУ заданному теоретическому закону."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Критерий согласия $\\chi^2$\n",
    "\n",
    "\n",
    "**Выборка**    $X_n = {X_1 ... X_n}$ из нормального распределения , при этом $\\sigma$ неизвестно.\n",
    "\n",
    "**Нулевая гипотеза**    $H_0: X \\sim N(\\mu,\\sigma^2)$\n",
    "\n",
    "**Альтернативная гипотеза**     $H_1: X \\not\\sim N(\\mu,\\sigma^2)$ \n",
    "\n",
    "Мы разбиваем область определения нашей величины на $K$ бинов (интервалов/бакетов/карманов). Границы этих интервалов – $a_i$. Подсчитаем следующие значения:\n",
    "\n",
    "$n_i$ - количество элементов нашей выборки в интервале $[a_i,a_{i+1}]$\n",
    "\n",
    "$p_i$ - вероятность попадания значения в интервал $[a_i,a_{i+1}]$ при условии справедливости $H_0$. Нам для этого надо будет взять функцию нормального распределения и посчитать ее разность в точках $a_i$ и $a_{i+1}$:\n",
    "\n",
    "$p_i = F_{N(\\mu,\\sigma^2)}(a_{i+1}) - F_{N(\\mu,\\sigma^2)}(a_{i})$\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "**Статистика** $\\chi^2(X_n) = \\sum^{K}_{i=1}\\dfrac{(n_i-np_i)^2}{np_i}$\n",
    "\n",
    "**Нулевое распределение** $\\chi^2(X_n) = \\begin{cases} \\chi^2_{K-1}, & \\text{Если нам известны $\\mu$ и $\\sigma$} \\\\ \\chi^2_{K-3}, & \\text{Если $\\mu$ и $\\sigma$ мы оцениваем по выборке} \\end{cases}$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Давайте для примера возьмем выборку из нормального распределения с параметрами $\\mu$ и $\\sigma$ и добавим к ней случайный шум $\\sim exponential(1)$. То есть наша случайная величина будет иметь распределение:\n",
    "$$X \\sim N(\\mu,\\sigma^2)+exponential(1)$$\n",
    "\n",
    "Пусть, например, "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "size = 5000\n",
    "bins = 20\n",
    "\n",
    "np.random.seed(31)\n",
    "sample = np.array([pylab.normal(70, 5) + pylab.exponential(1) for _ in range(size)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAigAAAGdCAYAAAA44ojeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAoDklEQVR4nO3df3DU9Z3H8Vd+LuHHbkwgu0QTflgVUhEVaLLqaU9zRBp7WuIPOimNJycjFzghipApgmJrGJw7Wu5E1HOAO0AtM0VLPMAQK15l5UcclKJEUGyCYRNOzC5gk0DyuT867LkCwoaE/WR5Pma+M+7389nd93vWr/vy+/3uJ3HGGCMAAACLxEe7AAAAgG8joAAAAOsQUAAAgHUIKAAAwDoEFAAAYB0CCgAAsA4BBQAAWIeAAgAArJMY7QI6o6OjQw0NDerXr5/i4uKiXQ4AADgHxhgdOXJEmZmZio//7nMkPTKgNDQ0KCsrK9plAACATqivr9dll132nXN6ZEDp16+fpL826HQ6o1wNAAA4F8FgUFlZWaHv8e/SIwPKycs6TqeTgAIAQA9zLrdncJMsAACwDgEFAABYh4ACAACsQ0ABAADWIaAAAADrEFAAAIB1CCgAAMA6BBQAAGCdiALK4MGDFRcXd8pWWloqSWppaVFpaanS09PVt29fFRUVqbGxMew16urqVFhYqN69eysjI0MzZ87UiRMnuq4jAADQ40UUULZv366DBw+GtqqqKknSPffcI0maMWOG1q1bpzVr1mjz5s1qaGjQ+PHjQ89vb29XYWGh2tratGXLFq1YsULLly/X3Llzu7AlAADQ08UZY0xnnzx9+nRVVlZq7969CgaDGjBggFavXq27775bkrRnzx4NHz5cPp9PeXl5Wr9+ve644w41NDTI7XZLkpYuXapZs2bp0KFDSk5OPqf3DQaDcrlcCgQCLHUPAEAPEcn3d6fvQWlra9PKlSv1wAMPKC4uTjU1NTp+/Ljy8/NDc4YNG6bs7Gz5fD5Jks/n04gRI0LhRJIKCgoUDAa1e/fuM75Xa2urgsFg2AYAAGJXpwPKa6+9pubmZt1///2SJL/fr+TkZKWmpobNc7vd8vv9oTnfDCcnx0+OnUlFRYVcLldoy8rK6mzZAACgB+h0QHnppZc0btw4ZWZmdmU9p1VeXq5AIBDa6uvru/09AQBA9CR25kl//vOftWnTJv3ud78L7fN4PGpra1Nzc3PYWZTGxkZ5PJ7QnG3btoW91slf+ZycczoOh0MOh6MzpQLogQbPfqPbXvvzBYXd9toAuk6nzqAsW7ZMGRkZKiz8/wN91KhRSkpKUnV1dWhfbW2t6urq5PV6JUler1e7du1SU1NTaE5VVZWcTqdycnI62wMAAIgxEZ9B6ejo0LJly1RSUqLExP9/usvl0qRJk1RWVqa0tDQ5nU5NmzZNXq9XeXl5kqSxY8cqJydHEydO1MKFC+X3+zVnzhyVlpZyhgQAAIREHFA2bdqkuro6PfDAA6eMLVq0SPHx8SoqKlJra6sKCgq0ZMmS0HhCQoIqKys1ZcoUeb1e9enTRyUlJZo/f/75dQEAAGLKea2DEi2sgwKcu+68n6Mn4h4UIHouyDooAAAA3YWAAgAArENAAQAA1iGgAAAA6xBQAACAdQgoAADAOgQUAABgHQIKAACwDgEFAABYh4ACAACsQ0ABAADWIaAAAADrEFAAAIB1CCgAAMA6BBQAAGAdAgoAALAOAQUAAFiHgAIAAKxDQAEAANYhoAAAAOsQUAAAgHUIKAAAwDoEFAAAYB0CCgAAsA4BBQAAWIeAAgAArENAAQAA1iGgAAAA6xBQAACAdQgoAADAOgQUAABgHQIKAACwDgEFAABYh4ACAACsQ0ABAADWIaAAAADrEFAAAIB1CCgAAMA6BBQAAGAdAgoAALBOxAHliy++0M9+9jOlp6crJSVFI0aM0I4dO0LjxhjNnTtXAwcOVEpKivLz87V3796w1zh8+LCKi4vldDqVmpqqSZMm6ejRo+ffDQAAiAkRBZSvvvpKN954o5KSkrR+/Xp99NFH+pd/+RddcskloTkLFy7U4sWLtXTpUm3dulV9+vRRQUGBWlpaQnOKi4u1e/duVVVVqbKyUu+8844mT57cdV0BAIAeLc4YY8518uzZs/Xuu+/qf/7nf047boxRZmamHnnkET366KOSpEAgILfbreXLl2vChAn6+OOPlZOTo+3bt2v06NGSpA0bNuhHP/qRDhw4oMzMzLPWEQwG5XK5FAgE5HQ6z7V84KI0ePYb0S7BKp8vKIx2CcBFK5Lv74jOoPz+97/X6NGjdc899ygjI0PXXXedXnzxxdD4/v375ff7lZ+fH9rncrmUm5srn88nSfL5fEpNTQ2FE0nKz89XfHy8tm7detr3bW1tVTAYDNsAAEDsiiigfPbZZ3ruued0xRVXaOPGjZoyZYr++Z//WStWrJAk+f1+SZLb7Q57ntvtDo35/X5lZGSEjScmJiotLS0059sqKirkcrlCW1ZWViRlAwCAHiaigNLR0aHrr79eTz/9tK677jpNnjxZDz74oJYuXdpd9UmSysvLFQgEQlt9fX23vh8AAIiuiALKwIEDlZOTE7Zv+PDhqqurkyR5PB5JUmNjY9icxsbG0JjH41FTU1PY+IkTJ3T48OHQnG9zOBxyOp1hGwAAiF0RBZQbb7xRtbW1Yfs++eQTDRo0SJI0ZMgQeTweVVdXh8aDwaC2bt0qr9crSfJ6vWpublZNTU1ozltvvaWOjg7l5uZ2uhEAABA7EiOZPGPGDN1www16+umnde+992rbtm164YUX9MILL0iS4uLiNH36dP3yl7/UFVdcoSFDhujxxx9XZmam7rrrLkl/PeNy++23hy4NHT9+XFOnTtWECRPO6Rc8AAAg9kUUUMaMGaO1a9eqvLxc8+fP15AhQ/TrX/9axcXFoTmPPfaYjh07psmTJ6u5uVk33XSTNmzYoF69eoXmrFq1SlOnTtVtt92m+Ph4FRUVafHixV3XFQAA6NEiWgfFFqyDApw71kEJxzooQPR02zooAAAAFwIBBQAAWIeAAgAArENAAQAA1iGgAAAA6xBQAACAdQgoAADAOgQUAABgHQIKAACwDgEFAABYh4ACAACsQ0ABAADWIaAAAADrEFAAAIB1CCgAAMA6BBQAAGAdAgoAALAOAQUAAFiHgAIAAKxDQAEAANYhoAAAAOsQUAAAgHUIKAAAwDoEFAAAYB0CCgAAsA4BBQAAWIeAAgAArENAAQAA1iGgAAAA6xBQAACAdQgoAADAOgQUAABgHQIKAACwDgEFAABYh4ACAACsQ0ABAADWIaAAAADrEFAAAIB1CCgAAMA6iZFMfuKJJ/Tkk0+G7bvqqqu0Z88eSVJLS4seeeQRvfLKK2ptbVVBQYGWLFkit9sdml9XV6cpU6boD3/4g/r27auSkhJVVFQoMTGiUoCYMnj2G9EuAQCsEnEq+P73v69Nmzb9/wt8I1jMmDFDb7zxhtasWSOXy6WpU6dq/PjxevfddyVJ7e3tKiwslMfj0ZYtW3Tw4EH9/Oc/V1JSkp5++ukuaAcAAMSCiANKYmKiPB7PKfsDgYBeeuklrV69WrfeeqskadmyZRo+fLjee+895eXl6c0339RHH32kTZs2ye1269prr9VTTz2lWbNm6YknnlBycvL5dwQAAHq8iO9B2bt3rzIzMzV06FAVFxerrq5OklRTU6Pjx48rPz8/NHfYsGHKzs6Wz+eTJPl8Po0YMSLskk9BQYGCwaB27959xvdsbW1VMBgM2wAAQOyKKKDk5uZq+fLl2rBhg5577jnt379ff/M3f6MjR47I7/crOTlZqampYc9xu93y+/2SJL/fHxZOTo6fHDuTiooKuVyu0JaVlRVJ2QAAoIeJ6BLPuHHjQv98zTXXKDc3V4MGDdJvf/tbpaSkdHlxJ5WXl6usrCz0OBgMElIAAIhh5/Uz49TUVF155ZXat2+fPB6P2tra1NzcHDansbExdM+Kx+NRY2PjKeMnx87E4XDI6XSGbQAAIHadV0A5evSoPv30Uw0cOFCjRo1SUlKSqqurQ+O1tbWqq6uT1+uVJHm9Xu3atUtNTU2hOVVVVXI6ncrJyTmfUgAAQAyJ6BLPo48+qh//+McaNGiQGhoaNG/ePCUkJOinP/2pXC6XJk2apLKyMqWlpcnpdGratGnyer3Ky8uTJI0dO1Y5OTmaOHGiFi5cKL/frzlz5qi0tFQOh6NbGgQAAD1PRAHlwIED+ulPf6ovv/xSAwYM0E033aT33ntPAwYMkCQtWrRI8fHxKioqCluo7aSEhARVVlZqypQp8nq96tOnj0pKSjR//vyu7QoAzqC7FsX7fEFht7wucLGKM8aYaBcRqWAwKJfLpUAgwP0oiAmsJNvzEVCAs4vk+5u/xQMAAKxDQAEAANYhoAAAAOsQUAAAgHUIKAAAwDoEFAAAYB0CCgAAsA4BBQAAWIeAAgAArENAAQAA1iGgAAAA6xBQAACAdQgoAADAOgQUAABgHQIKAACwDgEFAABYh4ACAACsQ0ABAADWIaAAAADrEFAAAIB1CCgAAMA6BBQAAGAdAgoAALAOAQUAAFiHgAIAAKxDQAEAANYhoAAAAOsQUAAAgHUIKAAAwDoEFAAAYB0CCgAAsA4BBQAAWCcx2gUAQCwYPPuNbnvtzxcUdttrA7biDAoAALAOAQUAAFiHgAIAAKxDQAEAANYhoAAAAOsQUAAAgHUIKAAAwDrnFVAWLFiguLg4TZ8+PbSvpaVFpaWlSk9PV9++fVVUVKTGxsaw59XV1amwsFC9e/dWRkaGZs6cqRMnTpxPKQAAIIZ0OqBs375dzz//vK655pqw/TNmzNC6deu0Zs0abd68WQ0NDRo/fnxovL29XYWFhWpra9OWLVu0YsUKLV++XHPnzu18FwAAIKZ0KqAcPXpUxcXFevHFF3XJJZeE9gcCAb300kv613/9V916660aNWqUli1bpi1btui9996TJL355pv66KOPtHLlSl177bUaN26cnnrqKT377LNqa2vrmq4AAECP1qmAUlpaqsLCQuXn54ftr6mp0fHjx8P2Dxs2TNnZ2fL5fJIkn8+nESNGyO12h+YUFBQoGAxq9+7dp32/1tZWBYPBsA0AAMSuiP8WzyuvvKL3339f27dvP2XM7/crOTlZqampYfvdbrf8fn9ozjfDycnxk2OnU1FRoSeffDLSUgEAQA8V0RmU+vp6Pfzww1q1apV69erVXTWdory8XIFAILTV19dfsPcGAAAXXkQBpaamRk1NTbr++uuVmJioxMREbd68WYsXL1ZiYqLcbrfa2trU3Nwc9rzGxkZ5PB5JksfjOeVXPScfn5zzbQ6HQ06nM2wDAACxK6KActttt2nXrl3auXNnaBs9erSKi4tD/5yUlKTq6urQc2pra1VXVyev1ytJ8nq92rVrl5qamkJzqqqq5HQ6lZOT00VtAQCAniyie1D69eunq6++Omxfnz59lJ6eHto/adIklZWVKS0tTU6nU9OmTZPX61VeXp4kaezYscrJydHEiRO1cOFC+f1+zZkzR6WlpXI4HF3UFgAA6Mkivkn2bBYtWqT4+HgVFRWptbVVBQUFWrJkSWg8ISFBlZWVmjJlirxer/r06aOSkhLNnz+/q0sBAAA9VJwxxkS7iEgFg0G5XC4FAgHuR0FMGDz7jWiXAIt9vqAw2iUAXSKS72/+Fg8AALAOAQUAAFiHgAIAAKxDQAEAANYhoAAAAOsQUAAAgHW6fB0UIJbxc2AAuDA4gwIAAKxDQAEAANYhoAAAAOsQUAAAgHUIKAAAwDoEFAAAYB0CCgAAsA4BBQAAWIeAAgAArENAAQAA1iGgAAAA6xBQAACAdQgoAADAOgQUAABgHQIKAACwDgEFAABYh4ACAACsQ0ABAADWIaAAAADrEFAAAIB1CCgAAMA6BBQAAGAdAgoAALAOAQUAAFiHgAIAAKxDQAEAANYhoAAAAOsQUAAAgHUIKAAAwDoEFAAAYB0CCgAAsA4BBQAAWCeigPLcc8/pmmuukdPplNPplNfr1fr160PjLS0tKi0tVXp6uvr27auioiI1NjaGvUZdXZ0KCwvVu3dvZWRkaObMmTpx4kTXdAMAAGJCRAHlsssu04IFC1RTU6MdO3bo1ltv1Z133qndu3dLkmbMmKF169ZpzZo12rx5sxoaGjR+/PjQ89vb21VYWKi2tjZt2bJFK1as0PLlyzV37tyu7QoAAPRoccYYcz4vkJaWpmeeeUZ33323BgwYoNWrV+vuu++WJO3Zs0fDhw+Xz+dTXl6e1q9frzvuuEMNDQ1yu92SpKVLl2rWrFk6dOiQkpOTz+k9g8GgXC6XAoGAnE7n+ZQPRGTw7DeiXQIuQp8vKIx2CUCXiOT7u9P3oLS3t+uVV17RsWPH5PV6VVNTo+PHjys/Pz80Z9iwYcrOzpbP55Mk+Xw+jRgxIhROJKmgoEDBYDB0FgYAACAx0ifs2rVLXq9XLS0t6tu3r9auXaucnBzt3LlTycnJSk1NDZvvdrvl9/slSX6/PyycnBw/OXYmra2tam1tDT0OBoORlg0AAHqQiM+gXHXVVdq5c6e2bt2qKVOmqKSkRB999FF31BZSUVEhl8sV2rKysrr1/QAAQHRFHFCSk5P1ve99T6NGjVJFRYVGjhyp3/zmN/J4PGpra1Nzc3PY/MbGRnk8HkmSx+M55Vc9Jx+fnHM65eXlCgQCoa2+vj7SsgEAQA9y3uugdHR0qLW1VaNGjVJSUpKqq6tDY7W1taqrq5PX65Ukeb1e7dq1S01NTaE5VVVVcjqdysnJOeN7OByO0E+bT24AACB2RXQPSnl5ucaNG6fs7GwdOXJEq1ev1ttvv62NGzfK5XJp0qRJKisrU1pampxOp6ZNmyav16u8vDxJ0tixY5WTk6OJEydq4cKF8vv9mjNnjkpLS+VwOLqlQQAA0PNEFFCampr085//XAcPHpTL5dI111yjjRs36u/+7u8kSYsWLVJ8fLyKiorU2tqqgoICLVmyJPT8hIQEVVZWasqUKfJ6verTp49KSko0f/78ru0KAAD0aOe9Dko0sA4KooV1UBANrIOCWHFB1kEBAADoLgQUAABgHQIKAACwDgEFAABYh4ACAACsQ0ABAADWIaAAAADrRPzXjAHbsVYJAPR8BBQAsFx3hW4WgIPNuMQDAACsQ0ABAADWIaAAAADrEFAAAIB1CCgAAMA6BBQAAGAdAgoAALAOAQUAAFiHgAIAAKxDQAEAANYhoAAAAOsQUAAAgHUIKAAAwDoEFAAAYB0CCgAAsA4BBQAAWIeAAgAArENAAQAA1iGgAAAA6xBQAACAdQgoAADAOgQUAABgHQIKAACwDgEFAABYh4ACAACsQ0ABAADWIaAAAADrEFAAAIB1CCgAAMA6BBQAAGAdAgoAALBORAGloqJCY8aMUb9+/ZSRkaG77rpLtbW1YXNaWlpUWlqq9PR09e3bV0VFRWpsbAybU1dXp8LCQvXu3VsZGRmaOXOmTpw4cf7dAACAmBBRQNm8ebNKS0v13nvvqaqqSsePH9fYsWN17Nix0JwZM2Zo3bp1WrNmjTZv3qyGhgaNHz8+NN7e3q7CwkK1tbVpy5YtWrFihZYvX665c+d2XVcAAKBHizPGmM4++dChQ8rIyNDmzZt18803KxAIaMCAAVq9erXuvvtuSdKePXs0fPhw+Xw+5eXlaf369brjjjvU0NAgt9stSVq6dKlmzZqlQ4cOKTk5+azvGwwG5XK5FAgE5HQ6O1s+YtTg2W9EuwSgR/h8QWG0S8BFJpLv7/O6ByUQCEiS0tLSJEk1NTU6fvy48vPzQ3OGDRum7Oxs+Xw+SZLP59OIESNC4USSCgoKFAwGtXv37tO+T2trq4LBYNgGAABiV6cDSkdHh6ZPn64bb7xRV199tSTJ7/crOTlZqampYXPdbrf8fn9ozjfDycnxk2OnU1FRIZfLFdqysrI6WzYAAOgBOh1QSktL9ac//UmvvPJKV9ZzWuXl5QoEAqGtvr6+298TAABET2JnnjR16lRVVlbqnXfe0WWXXRba7/F41NbWpubm5rCzKI2NjfJ4PKE527ZtC3u9k7/yOTnn2xwOhxwOR2dKBQAAPVBEZ1CMMZo6darWrl2rt956S0OGDAkbHzVqlJKSklRdXR3aV1tbq7q6Onm9XkmS1+vVrl271NTUFJpTVVUlp9OpnJyc8+kFAADEiIjOoJSWlmr16tV6/fXX1a9fv9A9Iy6XSykpKXK5XJo0aZLKysqUlpYmp9OpadOmyev1Ki8vT5I0duxY5eTkaOLEiVq4cKH8fr/mzJmj0tJSzpIAAABJEQaU5557TpL0wx/+MGz/smXLdP/990uSFi1apPj4eBUVFam1tVUFBQVasmRJaG5CQoIqKys1ZcoUeb1e9enTRyUlJZo/f/75dQIAAGLGea2DEi2sg4LvwjoowLlhHRRcaBdsHRQAAIDu0Klf8QAAer7uPNvI2RmcL86gAAAA6xBQAACAdQgoAADAOgQUAABgHQIKAACwDgEFAABYh4ACAACsQ0ABAADWIaAAAADrEFAAAIB1CCgAAMA6BBQAAGAdAgoAALAOAQUAAFiHgAIAAKxDQAEAANYhoAAAAOsQUAAAgHUIKAAAwDoEFAAAYB0CCgAAsA4BBQAAWIeAAgAArENAAQAA1iGgAAAA6xBQAACAdQgoAADAOgQUAABgHQIKAACwDgEFAABYh4ACAACsQ0ABAADWIaAAAADrEFAAAIB1CCgAAMA6BBQAAGAdAgoAALAOAQUAAFgn4oDyzjvv6Mc//rEyMzMVFxen1157LWzcGKO5c+dq4MCBSklJUX5+vvbu3Rs25/DhwyouLpbT6VRqaqomTZqko0ePnlcjAAAgdkQcUI4dO6aRI0fq2WefPe34woULtXjxYi1dulRbt25Vnz59VFBQoJaWltCc4uJi7d69W1VVVaqsrNQ777yjyZMnd74LAAAQU+KMMabTT46L09q1a3XXXXdJ+uvZk8zMTD3yyCN69NFHJUmBQEBut1vLly/XhAkT9PHHHysnJ0fbt2/X6NGjJUkbNmzQj370Ix04cECZmZlnfd9gMCiXy6VAICCn09nZ8hFlg2e/Ee0SAHSTzxcURrsEWCiS7+8uvQdl//798vv9ys/PD+1zuVzKzc2Vz+eTJPl8PqWmpobCiSTl5+crPj5eW7duPe3rtra2KhgMhm0AACB2dWlA8fv9kiS32x223+12h8b8fr8yMjLCxhMTE5WWlhaa820VFRVyuVyhLSsrqyvLBgAAlukRv+IpLy9XIBAIbfX19dEuCQAAdKMuDSgej0eS1NjYGLa/sbExNObxeNTU1BQ2fuLECR0+fDg059scDoecTmfYBgAAYleXBpQhQ4bI4/Gouro6tC8YDGrr1q3yer2SJK/Xq+bmZtXU1ITmvPXWW+ro6FBubm5XlgMAAHqoxEifcPToUe3bty/0eP/+/dq5c6fS0tKUnZ2t6dOn65e//KWuuOIKDRkyRI8//rgyMzNDv/QZPny4br/9dj344INaunSpjh8/rqlTp2rChAnn9AseAAAQ+yIOKDt27NDf/u3fhh6XlZVJkkpKSrR8+XI99thjOnbsmCZPnqzm5mbddNNN2rBhg3r16hV6zqpVqzR16lTddtttio+PV1FRkRYvXtwF7QAAgFhwXuugRAvroMQG1kEBYhfroOB0orYOCgAAQFeI+BIPAABn011nSDkzc/EgoOA7cRkGABANXOIBAADWIaAAAADrEFAAAIB1CCgAAMA6BBQAAGAdAgoAALAOAQUAAFiHgAIAAKxDQAEAANYhoAAAAOsQUAAAgHUIKAAAwDoEFAAAYB0CCgAAsA4BBQAAWIeAAgAArENAAQAA1kmMdgEAAJyrwbPf6LbX/nxBYbe9NiLHGRQAAGAdAgoAALAOAQUAAFiHgAIAAKxDQAEAANYhoAAAAOsQUAAAgHVYByVGdOfaAAAAXGicQQEAANYhoAAAAOsQUAAAgHW4BwUAAHXfvXz8jZ/O4QwKAACwDgEFAABYh4ACAACswz0oFxBrlQAAcG44gwIAAKxDQAEAANaJ6iWeZ599Vs8884z8fr9Gjhypf/u3f9MPfvCDaJYEAECX6s7L+7H8E+aonUF59dVXVVZWpnnz5un999/XyJEjVVBQoKampmiVBAAALBFnjDHReOPc3FyNGTNG//7v/y5J6ujoUFZWlqZNm6bZs2d/53ODwaBcLpcCgYCcTmeX18bNrACAnqCnnUGJ5Ps7Kpd42traVFNTo/Ly8tC++Ph45efny+fznTK/tbVVra2toceBQEDSXxvtDh2tX3fL6wIA0JWyZ6zpttf+05MFXf6aJ7+3z+XcSFQCyv/+7/+qvb1dbrc7bL/b7daePXtOmV9RUaEnn3zylP1ZWVndViMAABcz16+777WPHDkil8v1nXN6xDoo5eXlKisrCz3u6OjQ4cOHlZ6erri4uChWdqpgMKisrCzV19d3y+Unm12svV+sfUv0Tu/0fjHpit6NMTpy5IgyMzPPOjcqAaV///5KSEhQY2Nj2P7GxkZ5PJ5T5jscDjkcjrB9qamp3VnieXM6nRfdv7wnXay9X6x9S/RO7xcfeu9872c7c3JSVH7Fk5ycrFGjRqm6ujq0r6OjQ9XV1fJ6vdEoCQAAWCRql3jKyspUUlKi0aNH6wc/+IF+/etf69ixY/qHf/iHaJUEAAAsEbWAct999+nQoUOaO3eu/H6/rr32Wm3YsOGUG2d7GofDoXnz5p1ySepicLH2frH2LdE7vdP7xeRC9x61dVAAAADOhL/FAwAArENAAQAA1iGgAAAA6xBQAACAdQgonfDEE08oLi4ubBs2bFho/Ic//OEp4w899FAUK+5aX3zxhX72s58pPT1dKSkpGjFihHbs2BEaN8Zo7ty5GjhwoFJSUpSfn6+9e/dGseKuc7be77///lM++9tvvz2KFXeNwYMHn9JXXFycSktLJUktLS0qLS1Venq6+vbtq6KiolMWYuypztZ7rB7v7e3tevzxxzVkyBClpKTo8ssv11NPPRX2N1Ri9Vg/l95j9ViX/roM/fTp0zVo0CClpKTohhtu0Pbt20PjF+xzN4jYvHnzzPe//31z8ODB0Hbo0KHQ+C233GIefPDBsPFAIBDFirvO4cOHzaBBg8z9999vtm7daj777DOzceNGs2/fvtCcBQsWGJfLZV577TXzwQcfmL//+783Q4YMMX/5y1+iWPn5O5feS0pKzO233x722R8+fDiKVXeNpqamsJ6qqqqMJPOHP/zBGGPMQw89ZLKyskx1dbXZsWOHycvLMzfccEN0i+4iZ+s9Vo/3X/3qVyY9Pd1UVlaa/fv3mzVr1pi+ffua3/zmN6E5sXqsn0vvsXqsG2PMvffea3JycszmzZvN3r17zbx584zT6TQHDhwwxly4z52A0gnz5s0zI0eOPOP4LbfcYh5++OELVs+FNGvWLHPTTTedcbyjo8N4PB7zzDPPhPY1Nzcbh8NhXn755QtRYrc5W+/G/PU/WnfeeeeFKSiKHn74YXP55Zebjo4O09zcbJKSksyaNWtC4x9//LGRZHw+XxSr7B7f7N2Y2D3eCwsLzQMPPBC2b/z48aa4uNgYE9vH+tl6NyZ2j/Wvv/7aJCQkmMrKyrD9119/vfnFL35xQT93LvF00t69e5WZmamhQ4equLhYdXV1YeOrVq1S//79dfXVV6u8vFxff/11lCrtWr///e81evRo3XPPPcrIyNB1112nF198MTS+f/9++f1+5efnh/a5XC7l5ubK5/NFo+Quc7beT3r77beVkZGhq666SlOmTNGXX34ZhWq7T1tbm1auXKkHHnhAcXFxqqmp0fHjx8M+82HDhik7O7vHf+bf9u3eT4rF4/2GG25QdXW1PvnkE0nSBx98oD/+8Y8aN26cpNg+1s/W+0mxeKyfOHFC7e3t6tWrV9j+lJQU/fGPf7ywn3uXxp2LxH//93+b3/72t+aDDz4wGzZsMF6v12RnZ5tgMGiMMeb55583GzZsMB9++KFZuXKlufTSS81PfvKTKFfdNRwOh3E4HKa8vNy8//775vnnnze9evUyy5cvN8YY8+677xpJpqGhIex599xzj7n33nujUXKXOVvvxhjz8ssvm9dff918+OGHZu3atWb48OFmzJgx5sSJE1GsvGu9+uqrJiEhwXzxxRfGGGNWrVplkpOTT5k3ZswY89hjj13o8rrVt3s3JnaP9/b2djNr1iwTFxdnEhMTTVxcnHn66adD47F8rJ+td2Ni+1j3er3mlltuMV988YU5ceKE+a//+i8THx9vrrzyygv6uRNQusBXX31lnE6n+Y//+I/TjldXVxtJYfcq9FRJSUnG6/WG7Zs2bZrJy8szxsT2f7TO1vvpfPrpp0aS2bRpU3eXd8GMHTvW3HHHHaHHF1NA+XbvpxMrx/vLL79sLrvsMvPyyy+bDz/80Pznf/6nSUtLuyj+Z+RsvZ9OLB3r+/btMzfffLORZBISEsyYMWNMcXGxGTZs2AX93LnE0wVSU1N15ZVXat++facdz83NlaQzjvckAwcOVE5OTti+4cOHhy5xeTweSTrlFxyNjY2hsZ7qbL2fztChQ9W/f/+Y+Owl6c9//rM2bdqkf/zHfwzt83g8amtrU3Nzc9jcWPjMv+l0vZ9OrBzvM2fO1OzZszVhwgSNGDFCEydO1IwZM1RRUSEpto/1s/V+OrF0rF9++eXavHmzjh49qvr6em3btk3Hjx/X0KFDL+jnTkDpAkePHtWnn36qgQMHnnZ8586dknTG8Z7kxhtvVG1tbdi+Tz75RIMGDZIkDRkyRB6PR9XV1aHxYDCorVu3yuv1XtBau9rZej+dAwcO6Msvv4yJz16Sli1bpoyMDBUWFob2jRo1SklJSWGfeW1trerq6nr8Z/5Np+v9dGLleP/6668VHx/+FZGQkKCOjg5JsX2sn63304m1Y12S+vTpo4EDB+qrr77Sxo0bdeedd17Yz71Lz8dcJB555BHz9ttvm/3795t3333X5Ofnm/79+5umpiazb98+M3/+fLNjxw6zf/9+8/rrr5uhQ4eam2++Odpld4lt27aZxMRE86tf/crs3bvXrFq1yvTu3dusXLkyNGfBggUmNTU1dH32zjvvjImfHp6t9yNHjphHH33U+Hw+s3//frNp0yZz/fXXmyuuuMK0tLREufrz197ebrKzs82sWbNOGXvooYdMdna2eeutt8yOHTuM1+s95XJYT3am3mP5eC8pKTGXXnpp6Ke2v/vd70z//v3DLtvF6rF+tt5j/VjfsGGDWb9+vfnss8/Mm2++aUaOHGlyc3NNW1ubMebCfe4ElE647777zMCBA01ycrK59NJLzX333Re63lxXV2duvvlmk5aWZhwOh/ne975nZs6cGRPrIpy0bt06c/XVVxuHw2GGDRtmXnjhhbDxjo4O8/jjjxu3220cDoe57bbbTG1tbZSq7Vrf1fvXX39txo4dawYMGGCSkpLMoEGDzIMPPmj8fn8UK+46GzduNJJO+1n+5S9/Mf/0T/9kLrnkEtO7d2/zk5/8xBw8eDAKVXaPM/Uey8d7MBg0Dz/8sMnOzja9evUyQ4cONb/4xS9Ma2traE6sHutn6z3Wj/VXX33VDB061CQnJxuPx2NKS0tNc3NzaPxCfe5xxnxjaTwAAAALcA8KAACwDgEFAABYh4ACAACsQ0ABAADWIaAAAADrEFAAAIB1CCgAAMA6BBQAAGAdAgoAALAOAQUAAFiHgAIAAKxDQAEAANb5P3qH6bUi7BtTAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from matplotlib.pyplot import hist\n",
    "from scipy.stats import norm\n",
    "\n",
    "h1 = hist(sample, bins = bins)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "intervals = np.linspace(min(sample), max(sample), bins)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import chisquare"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "expected = []\n",
    "observed = []\n",
    "for i in range(1, len(intervals)):\n",
    "    observed.append(len(sample[(sample < intervals[i]) & (sample >= intervals[i - 1])]))\n",
    "    expected.append((norm.cdf(intervals[i], 70, 5) - norm.cdf(intervals[i - 1], 70, 5)) * len(sample))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4996.967652014366"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.array(expected).sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "expected = np.array(expected)\n",
    "expected[-1] = expected[-1] + size - expected.sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4999"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.array(observed).sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "observed = np.array(observed)\n",
    "observed[-1] = observed[-1] + size - observed.sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<BarContainer object of 19 artists>"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA9oAAAJGCAYAAACkxP3LAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAxh0lEQVR4nO3df3DU9Z348VcgEBRNUlASqYC0Z4u0qC1a2Np+WzFHpKnVkWnV4RSvnN4xwStytcoN/ij2xNKeejoR2w6CHWtt7VR7BYsCKt5JQIs6R6FDtUeFHibc1SNReyRAPt8/btjrKj9ceJNN9PGY2Rmyn/fuvj5+5pPM0092U5ZlWRYAAABAEn1KPQAAAAC8mwhtAAAASEhoAwAAQEJCGwAAABIS2gAAAJCQ0AYAAICEhDYAAAAkVF7qAQ5FV1dXbNu2LY499tgoKysr9TgAAAC8y2VZFq+//noMHTo0+vQ58DXrXhna27Zti2HDhpV6DAAAAN5jtm7dGieeeOIB1/TK0D722GMj4n93sLKyssTTAAAA8G7X3t4ew4YNy/fogfTK0N776+KVlZVCGwAAgG7zTt6+7MPQAAAAICGhDQAAAAkJbQAAAEhIaAMAAEBCQhsAAAASEtoAAACQkNAGAACAhIQ2AAAAJCS0AQAAICGhDQAAAAkJbQAAAEhIaAMAAEBCQhsAAAASEtoAAACQkNAGAACAhIQ2AAAAJCS0AQAAICGhDQAAAAkJbQAAAEhIaAMAAEBCQhsAAAASEtoAAACQkNAGAACAhMpLPQAA8A49Oa/UExQ6e3apJwCAHskVbQAAAEhIaAMAAEBCQhsAAAAS8h5tAN6TTrpuaalHKPC7WxtKPQIAkIgr2gAAAJCQ0AYAAICEhDYAAAAkJLQBAAAgIaENAAAACQltAAAASEhoAwAAQEJCGwAAABIS2gAAAJCQ0AYAAICEhDYAAAAkJLQBAAAgIaENAAAACQltAAAASEhoAwAAQEJCGwAAABIS2gAAAJCQ0AYAAICEhDYAAAAkVF7qAQCgFGaW/6TUI7xFQ6kHAAAScUUbAAAAEhLaAAAAkJDQBgAAgISENgAAACQktAEAACAhnzoOAL3EHSt/U+oRCsw8u9QTAEDP5Io2AAAAJCS0AQAAIKGiQvukk06KsrKyt90aGxsjImLnzp3R2NgYgwcPjmOOOSYmT54cra2tBc+xZcuWaGhoiKOPPjqGDBkS11xzTezevTvdHgEAAEAJFRXazz33XLz66qv52/LlyyMi4otf/GJERFx99dXx85//PB566KFYtWpVbNu2LS688ML84/fs2RMNDQ3R2dkZq1evjvvuuy8WL14cN9xwQ8JdAgAAgNIpKrSPP/74qK2tzd+WLFkSH/zgB+Mzn/lMtLW1xcKFC+O2226LCRMmxNixY2PRokWxevXqWLNmTUREPP7447Fx48a4//774/TTT49JkybFzTffHE1NTdHZ2XlEdhAAAAC60yG/R7uzszPuv//++PKXvxxlZWWxbt262LVrV9TV1eXXjBo1KoYPHx7Nzc0REdHc3BxjxoyJmpqa/Jr6+vpob2+PDRs27Pe1Ojo6or29veAGAAAAPdEhh/YjjzwSO3bsiMsvvzwiIlpaWqJ///5RXV1dsK6mpiZaWlrya/40svdu37ttf+bNmxdVVVX527Bhww51bAAAADiiDjm0Fy5cGJMmTYqhQ4emnGefZs+eHW1tbfnb1q1bj/hrAgAAwKEoP5QHvfLKK7FixYr46U9/mr+vtrY2Ojs7Y8eOHQVXtVtbW6O2tja/5tlnny14rr2fSr53zb5UVFRERUXFoYwKAAAA3eqQrmgvWrQohgwZEg0NDfn7xo4dG/369YuVK1fm79u0aVNs2bIlcrlcRETkcrlYv359bN++Pb9m+fLlUVlZGaNHjz7UfQAAAIAeo+gr2l1dXbFo0aKYOnVqlJf/38Orqqpi2rRpMWvWrBg0aFBUVlbGVVddFblcLsaPHx8RERMnTozRo0fHpZdeGvPnz4+WlpaYM2dONDY2umINAADAu0LRob1ixYrYsmVLfPnLX37btttvvz369OkTkydPjo6Ojqivr4+77747v71v376xZMmSmD59euRyuRg4cGBMnTo15s6de3h7AQAAAD1E0aE9ceLEyLJsn9sGDBgQTU1N0dTUtN/HjxgxIh599NFiXxYAAAB6hUP+1HEAAADg7YQ2AAAAJCS0AQAAIKFD+jvaAPCnTrpuaalHKPC7WxsOvggA4AhxRRsAAAASEtoAAACQkNAGAACAhIQ2AAAAJCS0AQAAICGhDQAAAAkJbQAAAEhIaAMAAEBCQhsAAAASEtoAAACQkNAGAACAhIQ2AAAAJCS0AQAAICGhDQAAAAkJbQAAAEhIaAMAAEBCQhsAAAASEtoAAACQkNAGAACAhIQ2AAAAJCS0AQAAICGhDQAAAAkJbQAAAEhIaAMAAEBCQhsAAAASEtoAAACQkNAGAACAhIQ2AAAAJCS0AQAAICGhDQAAAAkJbQAAAEhIaAMAAEBCQhsAAAASEtoAAACQkNAGAACAhIQ2AAAAJCS0AQAAICGhDQAAAAkJbQAAAEhIaAMAAEBCQhsAAAASEtoAAACQkNAGAACAhIQ2AAAAJCS0AQAAIKHyUg8AALyLPTmv1BMUOnt2qScA4D1AaANw2GaW/6TUI7xFQ6kHAADew/zqOAAAACQktAEAACAhoQ0AAAAJCW0AAABISGgDAABAQkIbAAAAEio6tP/jP/4j/uIv/iIGDx4cRx11VIwZMyZ++ctf5rdnWRY33HBDnHDCCXHUUUdFXV1dvPTSSwXP8dprr8WUKVOisrIyqqurY9q0afHGG28c/t4AAABAiRUV2v/93/8dZ511VvTr1y9+8YtfxMaNG+Mf//Ef433ve19+zfz58+POO++Me+65J9auXRsDBw6M+vr62LlzZ37NlClTYsOGDbF8+fJYsmRJPP3003HllVem2ysAAAAokfJiFn/zm9+MYcOGxaJFi/L3jRw5Mv/vLMvijjvuiDlz5sT5558fERHf//73o6amJh555JG4+OKL49e//nUsW7YsnnvuuTjjjDMiIuKuu+6Kz33uc/Htb387hg4dmmK/AAAAoCSKuqL9z//8z3HGGWfEF7/4xRgyZEh87GMfi+9973v57Zs3b46Wlpaoq6vL31dVVRXjxo2L5ubmiIhobm6O6urqfGRHRNTV1UWfPn1i7dq1+3zdjo6OaG9vL7gBAABAT1RUaP/7v/97LFiwIE4++eR47LHHYvr06fG3f/u3cd9990VEREtLS0RE1NTUFDyupqYmv62lpSWGDBlSsL28vDwGDRqUX/NW8+bNi6qqqvxt2LBhxYwNAAAA3aao0O7q6oqPf/zjccstt8THPvaxuPLKK+OKK66Ie+6550jNFxERs2fPjra2tvxt69atR/T1AAAA4FAVFdonnHBCjB49uuC+U045JbZs2RIREbW1tRER0draWrCmtbU1v622tja2b99esH337t3x2muv5de8VUVFRVRWVhbcAAAAoCcqKrTPOuus2LRpU8F9v/nNb2LEiBER8b8fjFZbWxsrV67Mb29vb4+1a9dGLpeLiIhcLhc7duyIdevW5dc88cQT0dXVFePGjTvkHQEAAICeoKhPHb/66qvjk5/8ZNxyyy3xpS99KZ599tn47ne/G9/97ncjIqKsrCxmzpwZ3/jGN+Lkk0+OkSNHxvXXXx9Dhw6NCy64ICL+9wr4ueeem/+V8127dsWMGTPi4osv9onjAAAA9HpFhfaZZ54ZDz/8cMyePTvmzp0bI0eOjDvuuCOmTJmSX/O1r30t3nzzzbjyyitjx44d8alPfSqWLVsWAwYMyK/5wQ9+EDNmzIhzzjkn+vTpE5MnT44777wz3V4BAABAiRQV2hERn//85+Pzn//8freXlZXF3LlzY+7cuftdM2jQoHjggQeKfWkAAADo8Yp6jzYAAABwYEIbAAAAEhLaAAAAkFDR79EGAHin7lj5m1KPUGDm2aWeAID3Ale0AQAAICGhDQAAAAkJbQAAAEhIaAMAAEBCQhsAAAASEtoAAACQkNAGAACAhIQ2AAAAJCS0AQAAICGhDQAAAAkJbQAAAEhIaAMAAEBCQhsAAAASEtoAAACQkNAGAACAhIQ2AAAAJCS0AQAAICGhDQAAAAkJbQAAAEhIaAMAAEBCQhsAAAASEtoAAACQkNAGAACAhIQ2AAAAJCS0AQAAICGhDQAAAAmVl3oAAAqddN3SUo9Q4He3NpR6BACAXsUVbQAAAEhIaAMAAEBCQhsAAAASEtoAAACQkNAGAACAhIQ2AAAAJCS0AQAAICGhDQAAAAkJbQAAAEhIaAMAAEBCQhsAAAASEtoAAACQkNAGAACAhIQ2AAAAJCS0AQAAICGhDQAAAAkJbQAAAEhIaAMAAEBCQhsAAAASEtoAAACQkNAGAACAhIQ2AAAAJCS0AQAAICGhDQAAAAkJbQAAAEhIaAMAAEBCRYX2TTfdFGVlZQW3UaNG5bfv3LkzGhsbY/DgwXHMMcfE5MmTo7W1teA5tmzZEg0NDXH00UfHkCFD4pprrondu3en2RsAAAAosfJiH/CRj3wkVqxY8X9PUP5/T3H11VfH0qVL46GHHoqqqqqYMWNGXHjhhfHMM89ERMSePXuioaEhamtrY/Xq1fHqq6/GZZddFv369Ytbbrklwe4AAABAaRUd2uXl5VFbW/u2+9va2mLhwoXxwAMPxIQJEyIiYtGiRXHKKafEmjVrYvz48fH444/Hxo0bY8WKFVFTUxOnn3563HzzzXHttdfGTTfdFP379z/8PQIAAIASKvo92i+99FIMHTo0PvCBD8SUKVNiy5YtERGxbt262LVrV9TV1eXXjho1KoYPHx7Nzc0REdHc3BxjxoyJmpqa/Jr6+vpob2+PDRs27Pc1Ozo6or29veAGAAAAPVFRoT1u3LhYvHhxLFu2LBYsWBCbN2+OT3/60/H6669HS0tL9O/fP6qrqwseU1NTEy0tLRER0dLSUhDZe7fv3bY/8+bNi6qqqvxt2LBhxYwNAAAA3aaoXx2fNGlS/t+nnnpqjBs3LkaMGBE//vGP46ijjko+3F6zZ8+OWbNm5b9ub28X2wAAAPRIh/Xnvaqrq+NDH/pQvPzyy1FbWxudnZ2xY8eOgjWtra3593TX1ta+7VPI9369r/d971VRURGVlZUFNwAAAOiJDiu033jjjfjtb38bJ5xwQowdOzb69esXK1euzG/ftGlTbNmyJXK5XERE5HK5WL9+fWzfvj2/Zvny5VFZWRmjR48+nFEAAACgRyjqV8e/+tWvxnnnnRcjRoyIbdu2xY033hh9+/aNSy65JKqqqmLatGkxa9asGDRoUFRWVsZVV10VuVwuxo8fHxEREydOjNGjR8ell14a8+fPj5aWlpgzZ040NjZGRUXFEdlBAAAA6E5Fhfbvf//7uOSSS+IPf/hDHH/88fGpT30q1qxZE8cff3xERNx+++3Rp0+fmDx5cnR0dER9fX3cfffd+cf37ds3lixZEtOnT49cLhcDBw6MqVOnxty5c9PuFQAAAJRIUaH94IMPHnD7gAEDoqmpKZqamva7ZsSIEfHoo48W87IAAADQaxzWe7QBAACAQkIbAAAAEhLaAAAAkJDQBgAAgISENgAAACQktAEAACAhoQ0AAAAJCW0AAABISGgDAABAQkIbAAAAEiov9QAAAD3Kk/NKPUGhs2eXegIAiuSKNgAAACQktAEAACAhoQ0AAAAJCW0AAABISGgDAABAQkIbAAAAEhLaAAAAkJDQBgAAgISENgAAACQktAEAACAhoQ0AAAAJCW0AAABISGgDAABAQkIbAAAAEhLaAAAAkJDQBgAAgISENgAAACRUXuoBACg0s/wnpR7hLRpKPQAAQK/iijYAAAAkJLQBAAAgIaENAAAACQltAAAASEhoAwAAQEJCGwAAABIS2gAAAJCQ0AYAAICEhDYAAAAkJLQBAAAgIaENAAAACQltAAAASEhoAwAAQEJCGwAAABIS2gAAAJCQ0AYAAICEhDYAAAAkJLQBAAAgIaENAAAACQltAAAASKi81AMAAPQkd6z8TalHKDDz7FJPAECxXNEGAACAhIQ2AAAAJCS0AQAAICGhDQAAAAkJbQAAAEhIaAMAAEBCQhsAAAASOqzQvvXWW6OsrCxmzpyZv2/nzp3R2NgYgwcPjmOOOSYmT54cra2tBY/bsmVLNDQ0xNFHHx1DhgyJa665Jnbv3n04owAAAECPcMih/dxzz8V3vvOdOPXUUwvuv/rqq+PnP/95PPTQQ7Fq1arYtm1bXHjhhfnte/bsiYaGhujs7IzVq1fHfffdF4sXL44bbrjh0PcCAAAAeohDCu033ngjpkyZEt/73vfife97X/7+tra2WLhwYdx2220xYcKEGDt2bCxatChWr14da9asiYiIxx9/PDZu3Bj3339/nH766TFp0qS4+eabo6mpKTo7O9PsFQAAAJTIIYV2Y2NjNDQ0RF1dXcH969ati127dhXcP2rUqBg+fHg0NzdHRERzc3OMGTMmampq8mvq6+ujvb09NmzYsM/X6+joiPb29oIbAAAA9ETlxT7gwQcfjOeffz6ee+65t21raWmJ/v37R3V1dcH9NTU10dLSkl/zp5G9d/vebfsyb968+PrXv17sqAAAANDtirqivXXr1vjKV74SP/jBD2LAgAFHaqa3mT17drS1teVvW7du7bbXBgAAgGIUFdrr1q2L7du3x8c//vEoLy+P8vLyWLVqVdx5551RXl4eNTU10dnZGTt27Ch4XGtra9TW1kZERG1t7ds+hXzv13vXvFVFRUVUVlYW3AAAAKAnKiq0zznnnFi/fn28+OKL+dsZZ5wRU6ZMyf+7X79+sXLlyvxjNm3aFFu2bIlcLhcREblcLtavXx/bt2/Pr1m+fHlUVlbG6NGjE+0WAAAAlEZR79E+9thj46Mf/WjBfQMHDozBgwfn7582bVrMmjUrBg0aFJWVlXHVVVdFLpeL8ePHR0TExIkTY/To0XHppZfG/Pnzo6WlJebMmRONjY1RUVGRaLcAAACgNIr+MLSDuf3226NPnz4xefLk6OjoiPr6+rj77rvz2/v27RtLliyJ6dOnRy6Xi4EDB8bUqVNj7ty5qUcBAACAbnfYof3UU08VfD1gwIBoamqKpqam/T5mxIgR8eijjx7uSwMAAECPc0h/RxsAAADYN6ENAAAACQltAAAASEhoAwAAQEJCGwAAABIS2gAAAJCQ0AYAAICEhDYAAAAkJLQBAAAgofJSDwBwJJ103dJSj1Dgd7c2lHoEAACOMFe0AQAAICGhDQAAAAkJbQAAAEhIaAMAAEBCQhsAAAASEtoAAACQkNAGAACAhIQ2AAAAJCS0AQAAICGhDQAAAAkJbQAAAEhIaAMAAEBCQhsAAAASEtoAAACQkNAGAACAhIQ2AAAAJCS0AQAAICGhDQAAAAkJbQAAAEhIaAMAAEBCQhsAAAASKi/1AAAAHKYn55V6gkJnzy71BAAl5Yo2AAAAJCS0AQAAICGhDQAAAAkJbQAAAEhIaAMAAEBCQhsAAAASEtoAAACQkNAGAACAhIQ2AAAAJCS0AQAAICGhDQAAAAkJbQAAAEhIaAMAAEBCQhsAAAASEtoAAACQkNAGAACAhIQ2AAAAJCS0AQAAICGhDQAAAAkJbQAAAEiovNQDABxJM8t/UuoR3qKh1AMAAHCEuaINAAAACQltAAAASEhoAwAAQEJCGwAAABIqKrQXLFgQp556alRWVkZlZWXkcrn4xS9+kd++c+fOaGxsjMGDB8cxxxwTkydPjtbW1oLn2LJlSzQ0NMTRRx8dQ4YMiWuuuSZ2796dZm8AAACgxIoK7RNPPDFuvfXWWLduXfzyl7+MCRMmxPnnnx8bNmyIiIirr746fv7zn8dDDz0Uq1atim3btsWFF16Yf/yePXuioaEhOjs7Y/Xq1XHffffF4sWL44Ybbki7VwAAAFAiRf15r/POO6/g63/4h3+IBQsWxJo1a+LEE0+MhQsXxgMPPBATJkyIiIhFixbFKaecEmvWrInx48fH448/Hhs3bowVK1ZETU1NnH766XHzzTfHtddeGzfddFP0798/3Z4BAABACRzye7T37NkTDz74YLz55puRy+Vi3bp1sWvXrqirq8uvGTVqVAwfPjyam5sjIqK5uTnGjBkTNTU1+TX19fXR3t6evyq+Lx0dHdHe3l5wAwAAgJ6o6NBev359HHPMMVFRURF/8zd/Ew8//HCMHj06Wlpaon///lFdXV2wvqamJlpaWiIioqWlpSCy927fu21/5s2bF1VVVfnbsGHDih0bAAAAukXRof3hD384XnzxxVi7dm1Mnz49pk6dGhs3bjwSs+XNnj072tra8retW7ce0dcDAACAQ1XUe7QjIvr37x9/9md/FhERY8eOjeeeey7+6Z/+KS666KLo7OyMHTt2FFzVbm1tjdra2oiIqK2tjWeffbbg+fZ+KvneNftSUVERFRUVxY4KAAAA3e6w/452V1dXdHR0xNixY6Nfv36xcuXK/LZNmzbFli1bIpfLRURELpeL9evXx/bt2/Nrli9fHpWVlTF69OjDHQUAAABKrqgr2rNnz45JkybF8OHD4/XXX48HHnggnnrqqXjssceiqqoqpk2bFrNmzYpBgwZFZWVlXHXVVZHL5WL8+PERETFx4sQYPXp0XHrppTF//vxoaWmJOXPmRGNjoyvWAAAAvCsUFdrbt2+Pyy67LF599dWoqqqKU089NR577LH48z//84iIuP3226NPnz4xefLk6OjoiPr6+rj77rvzj+/bt28sWbIkpk+fHrlcLgYOHBhTp06NuXPnpt0rAAAAKJGiQnvhwoUH3D5gwIBoamqKpqam/a4ZMWJEPProo8W8LAAAAPQah/0ebQAAAOD/CG0AAABISGgDAABAQkIbAAAAEhLaAAAAkJDQBgAAgISENgAAACQktAEAACAhoQ0AAAAJCW0AAABISGgDAABAQkIbAAAAEhLaAAAAkJDQBgAAgISENgAAACQktAEAACAhoQ0AAAAJCW0AAABISGgDAABAQkIbAAAAEiov9QAAAByeO1b+ptQjFJh5dqknACgtV7QBAAAgIaENAAAACQltAAAASEhoAwAAQEJCGwAAABIS2gAAAJCQ0AYAAICEhDYAAAAkJLQBAAAgIaENAAAACQltAAAASEhoAwAAQEJCGwAAABIS2gAAAJCQ0AYAAICEhDYAAAAkJLQBAAAgIaENAAAACQltAAAASEhoAwAAQEJCGwAAABIS2gAAAJCQ0AYAAICEyks9ANB7nHTd0lKPUOB3tzaUegQAAHgbV7QBAAAgIaENAAAACQltAAAASEhoAwAAQEJCGwAAABIS2gAAAJCQ0AYAAICEhDYAAAAkJLQBAAAgIaENAAAACQltAAAASEhoAwAAQEJFhfa8efPizDPPjGOPPTaGDBkSF1xwQWzatKlgzc6dO6OxsTEGDx4cxxxzTEyePDlaW1sL1mzZsiUaGhri6KOPjiFDhsQ111wTu3fvPvy9AQAAgBIrL2bxqlWrorGxMc4888zYvXt3/P3f/31MnDgxNm7cGAMHDoyIiKuvvjqWLl0aDz30UFRVVcWMGTPiwgsvjGeeeSYiIvbs2RMNDQ1RW1sbq1evjldffTUuu+yy6NevX9xyyy3p9xAAgJ7nyXmlnqDQ2bNLPQHwLlJUaC9btqzg68WLF8eQIUNi3bp18f/+3/+Ltra2WLhwYTzwwAMxYcKEiIhYtGhRnHLKKbFmzZoYP358PP7447Fx48ZYsWJF1NTUxOmnnx4333xzXHvttXHTTTdF//790+0dAAAAdLPDeo92W1tbREQMGjQoIiLWrVsXu3btirq6uvyaUaNGxfDhw6O5uTkiIpqbm2PMmDFRU1OTX1NfXx/t7e2xYcOGfb5OR0dHtLe3F9wAAACgJzrk0O7q6oqZM2fGWWedFR/96EcjIqKlpSX69+8f1dXVBWtramqipaUlv+ZPI3vv9r3b9mXevHlRVVWVvw0bNuxQxwYAAIAj6pBDu7GxMX71q1/Fgw8+mHKefZo9e3a0tbXlb1u3bj3irwkAAACHoqj3aO81Y8aMWLJkSTz99NNx4okn5u+vra2Nzs7O2LFjR8FV7dbW1qitrc2vefbZZwueb++nku9d81YVFRVRUVFxKKMCAABAtyrqinaWZTFjxox4+OGH44knnoiRI0cWbB87dmz069cvVq5cmb9v06ZNsWXLlsjlchERkcvlYv369bF9+/b8muXLl0dlZWWMHj36cPYFAAAASq6oK9qNjY3xwAMPxM9+9rM49thj8++prqqqiqOOOiqqqqpi2rRpMWvWrBg0aFBUVlbGVVddFblcLsaPHx8RERMnTozRo0fHpZdeGvPnz4+WlpaYM2dONDY2umoNAABAr1dUaC9YsCAiIj772c8W3L9o0aK4/PLLIyLi9ttvjz59+sTkyZOjo6Mj6uvr4+67786v7du3byxZsiSmT58euVwuBg4cGFOnTo25c+ce3p4AAABAD1BUaGdZdtA1AwYMiKampmhqatrvmhEjRsSjjz5azEsDAABAr3BYf0cbAAAAKCS0AQAAICGhDQAAAAkJbQAAAEhIaAMAAEBCQhsAAAASEtoAAACQkNAGAACAhIQ2AAAAJCS0AQAAICGhDQAAAAkJbQAAAEhIaAMAAEBCQhsAAAASEtoAAACQkNAGAACAhIQ2AAAAJCS0AQAAICGhDQAAAAkJbQAAAEhIaAMAAEBCQhsAAAASEtoAAACQUHmpBwB6j5nlPyn1CG/RUOoBAADgbVzRBgAAgISENgAAACQktAEAACAhoQ0AAAAJCW0AAABISGgDAABAQkIbAAAAEhLaAAAAkJDQBgAAgISENgAAACRUXuoBAAB477lj5W9KPUKBmWeXegLg3cQVbQAAAEhIaAMAAEBCQhsAAAASEtoAAACQkNAGAACAhIQ2AAAAJCS0AQAAICGhDQAAAAkJbQAAAEhIaAMAAEBCQhsAAAASEtoAAACQkNAGAACAhIQ2AAAAJFRe6gHgveqk65aWeoQCv7u1odQjAADAu4Ir2gAAAJCQ0AYAAICEhDYAAAAkJLQBAAAgIR+GBgAA78ST80o9QaGzZ5d6AmA/XNEGAACAhIQ2AAAAJFR0aD/99NNx3nnnxdChQ6OsrCweeeSRgu1ZlsUNN9wQJ5xwQhx11FFRV1cXL730UsGa1157LaZMmRKVlZVRXV0d06ZNizfeeOOwdgQAAAB6gqJD+80334zTTjstmpqa9rl9/vz5ceedd8Y999wTa9eujYEDB0Z9fX3s3Lkzv2bKlCmxYcOGWL58eSxZsiSefvrpuPLKKw99LwAAAKCHKPrD0CZNmhSTJk3a57Ysy+KOO+6IOXPmxPnnnx8REd///vejpqYmHnnkkbj44ovj17/+dSxbtiyee+65OOOMMyIi4q677orPfe5z8e1vfzuGDh16GLsDAAAApZX0PdqbN2+OlpaWqKury99XVVUV48aNi+bm5oiIaG5ujurq6nxkR0TU1dVFnz59Yu3atft83o6Ojmhvby+4AQAAQE+UNLRbWloiIqKmpqbg/pqamvy2lpaWGDJkSMH28vLyGDRoUH7NW82bNy+qqqryt2HDhqUcGwAAAJLpFZ86Pnv27Ghra8vftm7dWuqRAAAAYJ+ShnZtbW1ERLS2thbc39ramt9WW1sb27dvL9i+e/fueO211/Jr3qqioiIqKysLbgAAANATJQ3tkSNHRm1tbaxcuTJ/X3t7e6xduzZyuVxERORyudixY0esW7cuv+aJJ56Irq6uGDduXMpxAAAAoNsV/anjb7zxRrz88sv5rzdv3hwvvvhiDBo0KIYPHx4zZ86Mb3zjG3HyySfHyJEj4/rrr4+hQ4fGBRdcEBERp5xySpx77rlxxRVXxD333BO7du2KGTNmxMUXX+wTxwEAAOj1ig7tX/7yl3H22Wfnv541a1ZEREydOjUWL14cX/va1+LNN9+MK6+8Mnbs2BGf+tSnYtmyZTFgwID8Y37wgx/EjBkz4pxzzok+ffrE5MmT484770ywOwAAAFBaRYf2Zz/72ciybL/by8rKYu7cuTF37tz9rhk0aFA88MADxb40AAAA9Hi94lPHAQAAoLco+oo2kMbM8p+UeoS3aCj1AAAA8K7gijYAAAAkJLQBAAAgIaENAAAACQltAAAASEhoAwAAQEJCGwAAABIS2gAAAJCQ0AYAAICEhDYAAAAkJLQBAAAgIaENAAAACQltAAAASEhoAwAAQEJCGwAAABIS2gAAAJBQeakHAACA3uCOlb8p9QgFZp5d6gmA/XFFGwAAABIS2gAAAJCQ0AYAAICEhDYAAAAkJLQBAAAgIaENAAAACfnzXgAA8G715LxST1Do7NmlngC6hdDmXeGk65aWeoQCv7u1odQjAAAAJeJXxwEAACAhoQ0AAAAJCW0AAABISGgDAABAQkIbAAAAEhLaAAAAkJDQBgAAgISENgAAACQktAEAACAhoQ0AAAAJCW0AAABISGgDAABAQkIbAAAAEiov9QCQwszyn5R6hLdoKPUAAABAibiiDQAAAAkJbQAAAEjIr44DAAA9x5PzSj1BobNnl3oCeiGhDQAA71J3rPxNqUcoMPPsUk8A3cOvjgMAAEBCQhsAAAASEtoAAACQkPdo8zYnXbe01CMU+N2t/iY1AADQewht3mZm+U9KPcJbCG0AAKD38KvjAAAAkJDQBgAAgISENgAAACQktAEAACAhoQ0AAAAJCW0AAABIyJ/3AgAAeow7Vv6m1CMUmHn2O1j05LwjPkdRzp5d6gne84T2EXbSdUtLPUKB393qb1IDAAAcSSX91fGmpqY46aSTYsCAATFu3Lh49tlnSzkOAAAAHLaShfaPfvSjmDVrVtx4443x/PPPx2mnnRb19fWxffv2Uo0EAAAAh61kvzp+2223xRVXXBF/+Zd/GRER99xzTyxdujTuvffeuO6660o1VnIzy39S6hHewq+OAwBASr3xfeV3zPnLIz9IEWZ+Y1GpR0iqJKHd2dkZ69ati9mz/+9N+n369Im6urpobm5+2/qOjo7o6OjIf93W1hYREe3t7Ud+2MO0s6Oz1CMUeCf/zcx8+MzcPczcPczcPczcPczcPczcPczcPczcPXpD2+2dMcuyg64ty97JqsS2bdsW73//+2P16tWRy+Xy93/ta1+LVatWxdq1awvW33TTTfH1r3+9u8cEAACAAlu3bo0TTzzxgGt6xaeOz549O2bNmpX/uqurK1577bUYPHhwlJWVlXCy7tHe3h7Dhg2LrVu3RmVlZanH4RA4hr2b49f7OYa9m+PXuzl+vZ9j2Ls5fulkWRavv/56DB069KBrSxLaxx13XPTt2zdaW1sL7m9tbY3a2tq3ra+oqIiKioqC+6qrq4/kiD1SZWWlk6OXcwx7N8ev93MMezfHr3dz/Ho/x7B3c/zSqKqqekfrSvKp4/3794+xY8fGypUr8/d1dXXFypUrC36VHAAAAHqbkv3q+KxZs2Lq1KlxxhlnxCc+8Ym444474s0338x/CjkAAAD0RiUL7Ysuuij+8z//M2644YZoaWmJ008/PZYtWxY1NTWlGqnHqqioiBtvvPFtvz5P7+EY9m6OX+/nGPZujl/v5vj1fo5h7+b4lUZJPnUcAAAA3q1K8h5tAAAAeLcS2gAAAJCQ0AYAAICEhDYAAAAkJLQBAAAgIaHdQzQ1NcVJJ50UAwYMiHHjxsWzzz57wPUPPfRQjBo1KgYMGBBjxoyJRx99tJsm5a3mzZsXZ555Zhx77LExZMiQuOCCC2LTpk0HfMzixYujrKys4DZgwIBumpg/ddNNN73tWIwaNeqAj3H+9SwnnXTS245hWVlZNDY27nO986+0nn766TjvvPNi6NChUVZWFo888kjB9izL4oYbbogTTjghjjrqqKirq4uXXnrpoM9b7M9RDt2BjuGuXbvi2muvjTFjxsTAgQNj6NChcdlll8W2bdsO+JyH8r2YQ3Owc/Dyyy9/27E499xzD/q8zsHuc7BjuK+fiWVlZfGtb31rv8/pHExPaPcAP/rRj2LWrFlx4403xvPPPx+nnXZa1NfXx/bt2/e5fvXq1XHJJZfEtGnT4oUXXogLLrggLrjggvjVr37VzZMTEbFq1apobGyMNWvWxPLly2PXrl0xceLEePPNNw/4uMrKynj11Vfzt1deeaWbJuatPvKRjxQci3/913/d71rnX8/z3HPPFRy/5cuXR0TEF7/4xf0+xvlXOm+++Wacdtpp0dTUtM/t8+fPjzvvvDPuueeeWLt2bQwcODDq6+tj586d+33OYn+OcngOdAz/+Mc/xvPPPx/XX399PP/88/HTn/40Nm3aFF/4whcO+rzFfC/m0B3sHIyIOPfccwuOxQ9/+MMDPqdzsHsd7Bj+6bF79dVX4957742ysrKYPHnyAZ/XOZhYRsl94hOfyBobG/Nf79mzJxs6dGg2b968fa7/0pe+lDU0NBTcN27cuOyv//qvj+icvDPbt2/PIiJbtWrVftcsWrQoq6qq6r6h2K8bb7wxO+20097xeudfz/eVr3wl++AHP5h1dXXtc7vzr+eIiOzhhx/Of93V1ZXV1tZm3/rWt/L37dixI6uoqMh++MMf7vd5iv05SjpvPYb78uyzz2YRkb3yyiv7XVPs92LS2Nfxmzp1anb++ecX9TzOwdJ5J+fg+eefn02YMOGAa5yD6bmiXWKdnZ2xbt26qKury9/Xp0+fqKuri+bm5n0+prm5uWB9RER9ff1+19O92traIiJi0KBBB1z3xhtvxIgRI2LYsGFx/vnnx4YNG7pjPPbhpZdeiqFDh8YHPvCBmDJlSmzZsmW/a51/PVtnZ2fcf//98eUvfznKysr2u8751zNt3rw5WlpaCs6xqqqqGDdu3H7PsUP5OUr3amtri7Kysqiurj7gumK+F3NkPfXUUzFkyJD48Ic/HNOnT48//OEP+13rHOzZWltbY+nSpTFt2rSDrnUOpiW0S+y//uu/Ys+ePVFTU1Nwf01NTbS0tOzzMS0tLUWtp/t0dXXFzJkz46yzzoqPfvSj+1334Q9/OO6999742c9+Fvfff390dXXFJz/5yfj973/fjdMSETFu3LhYvHhxLFu2LBYsWBCbN2+OT3/60/H666/vc73zr2d75JFHYseOHXH55Zfvd43zr+faex4Vc44dys9Rus/OnTvj2muvjUsuuSQqKyv3u67Y78UcOeeee258//vfj5UrV8Y3v/nNWLVqVUyaNCn27Nmzz/XOwZ7tvvvui2OPPTYuvPDCA65zDqZXXuoB4N2ksbExfvWrXx30PS25XC5yuVz+609+8pNxyimnxHe+8524+eabj/SY/IlJkybl/33qqafGuHHjYsSIEfHjH//4Hf3fX3qWhQsXxqRJk2Lo0KH7XeP8g+6xa9eu+NKXvhRZlsWCBQsOuNb34p7j4osvzv97zJgxceqpp8YHP/jBeOqpp+Kcc84p4WQcinvvvTemTJly0A/9dA6m54p2iR133HHRt2/faG1tLbi/tbU1amtr9/mY2traotbTPWbMmBFLliyJJ598Mk488cSiHtuvX7/42Mc+Fi+//PIRmo53qrq6Oj70oQ/t91g4/3quV155JVasWBF/9Vd/VdTjnH89x97zqJhz7FB+jnLk7Y3sV155JZYvX37Aq9n7crDvxXSfD3zgA3Hcccft91g4B3uuf/mXf4lNmzYV/XMxwjmYgtAusf79+8fYsWNj5cqV+fu6urpi5cqVBVdc/lQulytYHxGxfPny/a7nyMqyLGbMmBEPP/xwPPHEEzFy5Miin2PPnj2xfv36OOGEE47AhBTjjTfeiN/+9rf7PRbOv55r0aJFMWTIkGhoaCjqcc6/nmPkyJFRW1tbcI61t7fH2rVr93uOHcrPUY6svZH90ksvxYoVK2Lw4MFFP8fBvhfTfX7/+9/HH/7wh/0eC+dgz7Vw4cIYO3ZsnHbaaUU/1jmYQKk/jY0se/DBB7OKiops8eLF2caNG7Mrr7wyq66uzlpaWrIsy7JLL700u+666/Lrn3nmmay8vDz79re/nf3617/Obrzxxqxfv37Z+vXrS7UL72nTp0/Pqqqqsqeeeip79dVX87c//vGP+TVvPYZf//rXs8ceeyz77W9/m61bty67+OKLswEDBmQbNmwoxS68p/3d3/1d9tRTT2WbN2/Onnnmmayuri477rjjsu3bt2dZ5vzrLfbs2ZMNHz48u/baa9+2zfnXs7z++uvZCy+8kL3wwgtZRGS33XZb9sILL+Q/kfrWW2/Nqqurs5/97GfZv/3bv2Xnn39+NnLkyOx//ud/8s8xYcKE7K677sp/fbCfo6R1oGPY2dmZfeELX8hOPPHE7MUXXyz4udjR0ZF/jrcew4N9LyadAx2/119/PfvqV7+aNTc3Z5s3b85WrFiRffzjH89OPvnkbOfOnfnncA6W1sG+j2ZZlrW1tWVHH310tmDBgn0+h3PwyBPaPcRdd92VDR8+POvfv3/2iU98IluzZk1+22c+85ls6tSpBet//OMfZx/60Iey/v37Zx/5yEeypUuXdvPE7BUR+7wtWrQov+atx3DmzJn5411TU5N97nOfy55//vnuH57soosuyk444YSsf//+2fvf//7soosuyl5++eX8dudf7/DYY49lEZFt2rTpbducfz3Lk08+uc/vmXuPUVdXV3b99ddnNTU1WUVFRXbOOee87biOGDEiu/HGGwvuO9DPUdI60DHcvHnzfn8uPvnkk/nneOsxPNj3YtI50PH74x//mE2cODE7/vjjs379+mUjRozIrrjiircFs3OwtA72fTTLsuw73/lOdtRRR2U7duzY53M4B4+8sizLsiN6yRwAAADeQ7xHGwAAABIS2gAAAJCQ0AYAAICEhDYAAAAkJLQBAAAgIaENAAAACQltAAAASEhoAwAAQEJCGwAAABIS2gAAAJCQ0AYAAICE/j8InZ4IhDYsEQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 1200x700 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "fig, ax = plt.subplots(1, 1, figsize = (12, 7))\n",
    "\n",
    "ax.bar(range(len(expected)), expected)\n",
    "ax.bar(range(len(observed)), observed, alpha = .5)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "253.73532675216688"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "((np.array(observed) - np.array(expected)) ** 2 / np.array(expected)).sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5000"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "observed.sum()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5000.0"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expected.sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Power_divergenceResult(statistic=253.73532675216688, pvalue=1.417119765427208e-43)"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "chisquare(observed, expected)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Наш тест уверенно отвергает нормальность с параметрами 70 и 5. В самом деле - распределение перекосило вправо. Давайте теперь оценим параметры нашего распределения по выборке. При этом кол-во степеней свободы уменьшается на 2. У распределения Хи-квадрат берется интеграл по правому хвосту от полученного в эксперименте значения статистики и до бесконечности."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "std = np.std(sample)\n",
    "mu = np.mean(sample)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(5.160590834308837, 71.02642515038143)"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "std, mu"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "expected = []\n",
    "observed = []\n",
    "for i in range(1, len(intervals)):\n",
    "    observed.append(len(sample[(sample < intervals[i]) & (sample >= intervals[i - 1])]))\n",
    "    expected.append((norm.cdf(intervals[i], mu, std) - norm.cdf(intervals[i-1], mu, std)) * len(sample))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "19"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(observed)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "14.171790988364513"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "((np.array(observed) - np.array(expected)) ** 2 / np.array(expected)).sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Power_divergenceResult(statistic=14.166096311096862, pvalue=0.5863414492872022)"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "chisquare(observed, np.array(expected) * np.mean(observed) / np.mean(expected), ddof = 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<BarContainer object of 19 artists>"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA9oAAAJGCAYAAACkxP3LAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAxY0lEQVR4nO3df3CU9Z3A8U8gEBRNUlQSqYC0Z1Va1BYtbG2vFXNEmlodmVYdTrHl9I4LXpH+UG7wF/bE0o5yOgHbDoIdS211qr2CRSEi3klAizpHoUe1R4UeJNzVI1F7JEj2/rhhr8tPF75kE329ZnaGPM93N5/1mSeZt8/upiSbzWYDAAAASKJXsQcAAACAdxOhDQAAAAkJbQAAAEhIaAMAAEBCQhsAAAASEtoAAACQkNAGAACAhEqLPcDh6OzsjK1bt8bxxx8fJSUlxR4HAACAd7lsNhtvvPFGDBo0KHr1Ovg16x4Z2lu3bo3BgwcXewwAAADeY7Zs2RKnnHLKQdf0yNA+/vjjI+L/nmB5eXmRpwEAAODdrq2tLQYPHpzr0YPpkaG95+Xi5eXlQhsAAIAu807evuzD0AAAACAhoQ0AAAAJCW0AAABISGgDAABAQkIbAAAAEhLaAAAAkJDQBgAAgISENgAAACQktAEAACAhoQ0AAAAJCW0AAABISGgDAABAQkIbAAAAEhLaAAAAkJDQBgAAgISENgAAACQktAEAACAhoQ0AAAAJCW0AAABISGgDAABAQkIbAAAAEhLaAAAAkJDQBgAAgIRKiz0AAPAOrZhV7AnyXTC92BMAQLfkijYAAAAkJLQBAAAgIaENAAAACQltAAAASEhoAwAAQEJCGwAAABIS2gAAAJCQ0AYAAICEhDYAAAAkJLQBAAAgIaENAAAACQltAAAASEhoAwAAQEJCGwAAABIS2gAAAJCQ0AYAAICEhDYAAAAkJLQBAAAgIaENAAAACQltAAAASEhoAwAAQEJCGwAAABIqLfYAAMA7M6fxN8UeIc/UC4o9AQB0T0IbgPemFbOKPUG+C6YXewIAIBEvHQcAAICEhDYAAAAkJLQBAAAgIaENAAAACQltAAAASEhoAwAAQEJCGwAAABIS2gAAAJBQaSGLTz311Hjttdf22f63f/u30dDQEDt37oyvfvWr8fDDD0d7e3vU1tbG3Llzo6qqKrd28+bNMXny5FixYkUcd9xxMXHixJg1a1aUlhY0CgAckTmNvyn2CHmmXlDsCQCAVAq6ov3CCy/Etm3bcrdly5ZFRMQXvvCFiIi44YYb4uc//3k88sgjsXLlyti6dWtcdtllufvv3r076urqoqOjI1atWhUPPvhgLFy4MG655ZaETwkAAACKp6DQPumkk6K6ujp3W7x4cXzwgx+MT3/609Ha2hrz58+Pu+++O8aMGRMjR46MBQsWxKpVq2L16tUREfHUU0/Fhg0b4qGHHopzzjknxo0bF3fccUc0NDRER0fHUXmCAAAA0JUO+z3aHR0d8dBDD8WXv/zlKCkpibVr18auXbuipqYmt+aMM86IIUOGRFNTU0RENDU1xYgRI/JeSl5bWxttbW2xfv36A36v9vb2aGtry7sBAABAd3TYof3444/Hjh074pprromIiObm5ujbt29UVlbmrauqqorm5ubcmj+N7D379+w7kFmzZkVFRUXuNnjw4MMdGwAAAI6qww7t+fPnx7hx42LQoEEp59mv6dOnR2tra+62ZcuWo/49AQAA4HAc1kd9v/baa7F8+fL46U9/mttWXV0dHR0dsWPHjryr2i0tLVFdXZ1b8/zzz+c9VktLS27fgZSVlUVZWdnhjAoAAABd6rCuaC9YsCAGDhwYdXV1uW0jR46MPn36RGNjY27bxo0bY/PmzZHJZCIiIpPJxLp162L79u25NcuWLYvy8vIYPnz44T4HAAAA6DYKvqLd2dkZCxYsiIkTJ+b97euKioqYNGlSTJs2LQYMGBDl5eVx/fXXRyaTidGjR0dExNixY2P48OFx1VVXxezZs6O5uTlmzJgR9fX1rlgDAADwrlBwaC9fvjw2b94cX/7yl/fZd88990SvXr1i/Pjx0d7eHrW1tTF37tzc/t69e8fixYtj8uTJkclkon///jFx4sSYOXPmkT0LAAAA6CYKDu2xY8dGNpvd775+/fpFQ0NDNDQ0HPD+Q4cOjSeeeKLQbwsAAAA9wmF/6jgAAACwL6ENAAAACQltAAAASEhoAwAAQEJCGwAAABIS2gAAAJCQ0AYAAICEhDYAAAAkJLQBAAAgIaENAAAACQltAAAASEhoAwAAQEJCGwAAABIS2gAAAJCQ0AYAAICEhDYAAAAkJLQBAAAgIaENAAAACQltAAAASEhoAwAAQEJCGwAAABIS2gAAAJCQ0AYAAICEhDYAAAAkJLQBAAAgIaENAAAACQltAAAASEhoAwAAQEJCGwAAABIS2gAAAJBQabEHAADevU69aUmxR8jzu7vqij0CAO8BrmgDAABAQkIbAAAAEhLaAAAAkJDQBgAAgISENgAAACQktAEAACAhoQ0AAAAJCW0AAABISGgDAABAQkIbAAAAEhLaAAAAkJDQBgAAgISENgAAACQktAEAACAhoQ0AAAAJCW0AAABISGgDAABAQkIbAAAAEiot9gAAwLvX1NJHiz3CXuqKPQAA7wFCG4Ajt2JWsSfId8H0Yk8AALyHeek4AAAAJCS0AQAAICGhDQAAAAkJbQAAAEhIaAMAAEBCQhsAAAASKji0/+M//iP+8i//Mk444YQ45phjYsSIEfHLX/4ytz+bzcYtt9wSJ598chxzzDFRU1MTr7zySt5jvP766zFhwoQoLy+PysrKmDRpUrz55ptH/mwAAACgyAr6O9r//d//Heeff35ccMEF8Ytf/CJOOumkeOWVV+J973tfbs3s2bPj3nvvjQcffDCGDRsWN998c9TW1saGDRuiX79+ERExYcKE2LZtWyxbtix27doVX/rSl+K6666LRYsWpX12AHSJOY2/KfYIeaZeUOwJAID3soJC+1vf+lYMHjw4FixYkNs2bNiw3L+z2WzMmTMnZsyYEZdccklERPzgBz+IqqqqePzxx+OKK66IX//617F06dJ44YUX4txzz42IiPvuuy8++9nPxne+850YNGhQiucFAAAARVHQS8f/6Z/+Kc4999z4whe+EAMHDoyPfvSj8f3vfz+3f9OmTdHc3Bw1NTW5bRUVFTFq1KhoamqKiIimpqaorKzMRXZERE1NTfTq1SvWrFmz3+/b3t4ebW1teTcAAADojgoK7X//93+PefPmxWmnnRZPPvlkTJ48Of7u7/4uHnzwwYiIaG5ujoiIqqqqvPtVVVXl9jU3N8fAgQPz9peWlsaAAQNya/Y2a9asqKioyN0GDx5cyNgAAADQZQoK7c7OzvjYxz4Wd955Z3z0ox+N6667Lq699tq4//77j9Z8ERExffr0aG1tzd22bNlyVL8fAAAAHK6CQvvkk0+O4cOH520788wzY/PmzRERUV1dHRERLS0teWtaWlpy+6qrq2P79u15+99+++14/fXXc2v2VlZWFuXl5Xk3AAAA6I4KCu3zzz8/Nm7cmLftN7/5TQwdOjQi/u+D0aqrq6OxsTG3v62tLdasWROZTCYiIjKZTOzYsSPWrl2bW/P0009HZ2dnjBo16rCfCAAAAHQHBX3q+A033BCf+MQn4s4774wvfvGL8fzzz8f3vve9+N73vhcRESUlJTF16tT45je/Gaeddlruz3sNGjQoLr300oj4vyvgF110Ue4l57t27YopU6bEFVdc4RPHAQAA6PEKCu3zzjsvHnvssZg+fXrMnDkzhg0bFnPmzIkJEybk1nzjG9+It956K6677rrYsWNHfPKTn4ylS5fm/oZ2RMQPf/jDmDJlSlx44YXRq1evGD9+fNx7773pnhUAAAAUSUGhHRHxuc99Lj73uc8dcH9JSUnMnDkzZs6cecA1AwYMiEWLFhX6rQEAAKDbK+g92gAAAMDBCW0AAABISGgDAABAQkIbAAAAEhLaAAAAkJDQBgAAgISENgAAACQktAEAACAhoQ0AAAAJCW0AAABISGgDAABAQkIbAAAAEhLaAAAAkJDQBgAAgISENgAAACQktAEAACAhoQ0AAAAJCW0AAABISGgDAABAQkIbAAAAEhLaAAAAkJDQBgAAgISENgAAACQktAEAACAhoQ0AAAAJCW0AAABISGgDAABAQkIbAAAAEhLaAAAAkJDQBgAAgISENgAAACQktAEAACAhoQ0AAAAJCW0AAABISGgDAABAQkIbAAAAEhLaAAAAkJDQBgAAgISENgAAACQktAEAACAhoQ0AAAAJCW0AAABISGgDAABAQkIbAAAAEhLaAAAAkJDQBgAAgISENgAAACQktAEAACAhoQ0AAAAJCW0AAABISGgDAABAQkIbAAAAEhLaAAAAkJDQBgAAgISENgAAACQktAEAACAhoQ0AAAAJlRay+Lbbbovbb789b9vpp58e//Zv/xYRETt37oyvfvWr8fDDD0d7e3vU1tbG3Llzo6qqKrd+8+bNMXny5FixYkUcd9xxMXHixJg1a1aUlhY0CgDAUXHqTUuKPUKe391VV+wRAChQwXX74Q9/OJYvX/7/D/AngXzDDTfEkiVL4pFHHomKioqYMmVKXHbZZfHcc89FRMTu3bujrq4uqqurY9WqVbFt27a4+uqro0+fPnHnnXcmeDoAAABQXAWHdmlpaVRXV++zvbW1NebPnx+LFi2KMWPGRETEggUL4swzz4zVq1fH6NGj46mnnooNGzbE8uXLo6qqKs4555y444474sYbb4zbbrst+vbte+TPCAAAAIqo4Pdov/LKKzFo0KD4wAc+EBMmTIjNmzdHRMTatWtj165dUVNTk1t7xhlnxJAhQ6KpqSkiIpqammLEiBF5LyWvra2Ntra2WL9+/QG/Z3t7e7S1teXdAAAAoDsqKLRHjRoVCxcujKVLl8a8efNi06ZN8alPfSreeOONaG5ujr59+0ZlZWXefaqqqqK5uTkiIpqbm/Mie8/+PfsOZNasWVFRUZG7DR48uJCxAQAAoMsU9NLxcePG5f591llnxahRo2Lo0KHxk5/8JI455pjkw+0xffr0mDZtWu7rtrY2sQ0AAEC3dER/3quysjI+9KEPxauvvhrV1dXR0dERO3bsyFvT0tKSe093dXV1tLS07LN/z74DKSsri/Ly8rwbAAAAdEdHFNpvvvlm/Pa3v42TTz45Ro4cGX369InGxsbc/o0bN8bmzZsjk8lEREQmk4l169bF9u3bc2uWLVsW5eXlMXz48CMZBQAAALqFgl46/rWvfS0uvvjiGDp0aGzdujVuvfXW6N27d1x55ZVRUVERkyZNimnTpsWAAQOivLw8rr/++shkMjF69OiIiBg7dmwMHz48rrrqqpg9e3Y0NzfHjBkzor6+PsrKyo7KEwQAAICuVFBo//73v48rr7wy/vCHP8RJJ50Un/zkJ2P16tVx0kknRUTEPffcE7169Yrx48dHe3t71NbWxty5c3P37927dyxevDgmT54cmUwm+vfvHxMnToyZM2emfVYAAABQJAWF9sMPP3zQ/f369YuGhoZoaGg44JqhQ4fGE088Uci3BQAAgB7jiN6jDQAAAOQT2gAAAJCQ0AYAAICEhDYAAAAkJLQBAAAgIaENAAAACQltAAAASEhoAwAAQEJCGwAAABIS2gAAAJBQabEHAADoTqaWPlrsEfZSV+wBACiQK9oAAACQkNAGAACAhIQ2AAAAJCS0AQAAICGhDQAAAAkJbQAAAEhIaAMAAEBCQhsAAAASEtoAAACQkNAGAACAhIQ2AAAAJCS0AQAAICGhDQAAAAkJbQAAAEiotNgDAJDv1JuWFHuEPL+7q67YIwAA9CiuaAMAAEBCQhsAAAAS8tJxgG5maumjxR5hL146DgBQCFe0AQAAICGhDQAAAAkJbQAAAEhIaAMAAEBCQhsAAAASEtoAAACQkNAGAACAhIQ2AAAAJCS0AQAAICGhDQAAAAkJbQAAAEhIaAMAAEBCQhsAAAASEtoAAACQkNAGAACAhIQ2AAAAJCS0AQAAICGhDQAAAAkJbQAAAEhIaAMAAEBCQhsAAAASEtoAAACQkNAGAACAhIQ2AAAAJCS0AQAAICGhDQAAAAkJbQAAAEhIaAMAAEBCRxTad911V5SUlMTUqVNz23bu3Bn19fVxwgknxHHHHRfjx4+PlpaWvPtt3rw56urq4thjj42BAwfG17/+9Xj77bePZBQAAADoFg47tF944YX47ne/G2eddVbe9htuuCF+/vOfxyOPPBIrV66MrVu3xmWXXZbbv3v37qirq4uOjo5YtWpVPPjgg7Fw4cK45ZZbDv9ZAAAAQDdxWKH95ptvxoQJE+L73/9+vO9978ttb21tjfnz58fdd98dY8aMiZEjR8aCBQti1apVsXr16oiIeOqpp2LDhg3x0EMPxTnnnBPjxo2LO+64IxoaGqKjoyPNswIAAIAiOazQrq+vj7q6uqipqcnbvnbt2ti1a1fe9jPOOCOGDBkSTU1NERHR1NQUI0aMiKqqqtya2traaGtri/Xr1+/3+7W3t0dbW1veDQAAALqj0kLv8PDDD8eLL74YL7zwwj77mpubo2/fvlFZWZm3vaqqKpqbm3Nr/jSy9+zfs29/Zs2aFbfffnuhowIAAECXK+iK9pYtW+IrX/lK/PCHP4x+/fodrZn2MX369Ghtbc3dtmzZ0mXfGwAAAApRUGivXbs2tm/fHh/72MeitLQ0SktLY+XKlXHvvfdGaWlpVFVVRUdHR+zYsSPvfi0tLVFdXR0REdXV1ft8Cvmer/es2VtZWVmUl5fn3QAAAKA7Kii0L7zwwli3bl28/PLLudu5554bEyZMyP27T58+0djYmLvPxo0bY/PmzZHJZCIiIpPJxLp162L79u25NcuWLYvy8vIYPnx4oqcFAAAAxVHQe7SPP/74+MhHPpK3rX///nHCCSfktk+aNCmmTZsWAwYMiPLy8rj++usjk8nE6NGjIyJi7NixMXz48Ljqqqti9uzZ0dzcHDNmzIj6+vooKytL9LQAAACgOAr+MLRDueeee6JXr14xfvz4aG9vj9ra2pg7d25uf+/evWPx4sUxefLkyGQy0b9//5g4cWLMnDkz9SgAAADQ5Y44tJ955pm8r/v16xcNDQ3R0NBwwPsMHTo0nnjiiSP91gAAANDtHNbf0QYAAAD2T2gDAABAQkIbAAAAEhLaAAAAkJDQBgAAgISENgAAACQktAEAACAhoQ0AAAAJCW0AAABISGgDAABAQkIbAAAAEhLaAAAAkJDQBgAAgISENgAAACQktAEAACAhoQ0AAAAJCW0AAABISGgDAABAQkIbAAAAEhLaAAAAkJDQBgAAgISENgAAACQktAEAACAhoQ0AAAAJCW0AAABISGgDAABAQkIbAAAAEhLaAAAAkFBpsQcAAODInHrTkmKPkOd3d9UVewSAonJFGwAAABIS2gAAAJCQl44DAPRwU0sfLfYIe/HSceC9zRVtAAAASEhoAwAAQEJCGwAAABIS2gAAAJCQ0AYAAICEhDYAAAAkJLQBAAAgIaENAAAACQltAAAASEhoAwAAQEJCGwAAABIS2gAAAJCQ0AYAAICEhDYAAAAkJLQBAAAgIaENAAAACQltAAAASEhoAwAAQEJCGwAAABIS2gAAAJBQabEHADiqVswq9gT5Lphe7AkAADjKXNEGAACAhIQ2AAAAJCS0AQAAICGhDQAAAAkVFNrz5s2Ls846K8rLy6O8vDwymUz84he/yO3fuXNn1NfXxwknnBDHHXdcjB8/PlpaWvIeY/PmzVFXVxfHHntsDBw4ML7+9a/H22+/nebZAAAAQJEVFNqnnHJK3HXXXbF27dr45S9/GWPGjIlLLrkk1q9fHxERN9xwQ/z85z+PRx55JFauXBlbt26Nyy67LHf/3bt3R11dXXR0dMSqVaviwQcfjIULF8Ytt9yS9lkBAABAkRT0570uvvjivK//4R/+IebNmxerV6+OU045JebPnx+LFi2KMWPGRETEggUL4swzz4zVq1fH6NGj46mnnooNGzbE8uXLo6qqKs4555y444474sYbb4zbbrst+vbtm+6ZAQAAQBEc9nu0d+/eHQ8//HC89dZbkclkYu3atbFr166oqanJrTnjjDNiyJAh0dTUFBERTU1NMWLEiKiqqsqtqa2tjba2ttxV8f1pb2+Ptra2vBsAAAB0RwWH9rp16+K4446LsrKy+Ju/+Zt47LHHYvjw4dHc3Bx9+/aNysrKvPVVVVXR3NwcERHNzc15kb1n/559BzJr1qyoqKjI3QYPHlzo2AAAANAlCnrpeETE6aefHi+//HK0trbGo48+GhMnToyVK1cejdlypk+fHtOmTct93dbWJraBd2RO42+KPUKeqRcUewIAAI62gkO7b9++8Wd/9mcRETFy5Mh44YUX4h//8R/j8ssvj46OjtixY0feVe2Wlpaorq6OiIjq6up4/vnn8x5vz6eS71mzP2VlZVFWVlboqAAAANDljvjvaHd2dkZ7e3uMHDky+vTpE42Njbl9GzdujM2bN0cmk4mIiEwmE+vWrYvt27fn1ixbtizKy8tj+PDhRzoKAAAAFF1BV7SnT58e48aNiyFDhsQbb7wRixYtimeeeSaefPLJqKioiEmTJsW0adNiwIABUV5eHtdff31kMpkYPXp0RESMHTs2hg8fHldddVXMnj07mpubY8aMGVFfX++KNQAAAO8KBYX29u3b4+qrr45t27ZFRUVFnHXWWfHkk0/GX/zFX0RExD333BO9evWK8ePHR3t7e9TW1sbcuXNz9+/du3csXrw4Jk+eHJlMJvr37x8TJ06MmTNnpn1WAAAAUCQFhfb8+fMPur9fv37R0NAQDQ0NB1wzdOjQeOKJJwr5tgAAANBjHPF7tAEAAID/J7QBAAAgIaENAAAACQltAAAASEhoAwAAQEJCGwAAABIS2gAAAJCQ0AYAAICEhDYAAAAkJLQBAAAgIaENAAAACQltAAAASEhoAwAAQEJCGwAAABIS2gAAAJCQ0AYAAICEhDYAAAAkJLQBAAAgIaENAAAACQltAAAASEhoAwAAQEJCGwAAABIS2gAAAJCQ0AYAAICEhDYAAAAkJLQBAAAgIaENAAAACQltAAAASEhoAwAAQEJCGwAAABIS2gAAAJCQ0AYAAICEhDYAAAAkJLQBAAAgIaENAAAACQltAAAASEhoAwAAQEJCGwAAABIS2gAAAJCQ0AYAAICEhDYAAAAkJLQBAAAgIaENAAAACQltAAAASEhoAwAAQEJCGwAAABIS2gAAAJCQ0AYAAICEhDYAAAAkJLQBAAAgIaENAAAACQltAAAASEhoAwAAQEJCGwAAABIS2gAAAJBQabEHAADgvefUm5YUe4Q8v7urrtgjAO8irmgDAABAQgWF9qxZs+K8886L448/PgYOHBiXXnppbNy4MW/Nzp07o76+Pk444YQ47rjjYvz48dHS0pK3ZvPmzVFXVxfHHntsDBw4ML7+9a/H22+/feTPBgAAAIqsoJeOr1y5Murr6+O8886Lt99+O/7+7/8+xo4dGxs2bIj+/ftHRMQNN9wQS5YsiUceeSQqKipiypQpcdlll8Vzzz0XERG7d++Ourq6qK6ujlWrVsW2bdvi6quvjj59+sSdd96Z/hkCANDtTC19tNgj7MVLx4F0CgrtpUuX5n29cOHCGDhwYKxduzb+/M//PFpbW2P+/PmxaNGiGDNmTERELFiwIM4888xYvXp1jB49Op566qnYsGFDLF++PKqqquKcc86JO+64I2688ca47bbbom/fvumeHQAAAHSxI3qPdmtra0REDBgwICIi1q5dG7t27YqamprcmjPOOCOGDBkSTU1NERHR1NQUI0aMiKqqqtya2traaGtri/Xr1+/3+7S3t0dbW1veDQAAALqjww7tzs7OmDp1apx//vnxkY98JCIimpubo2/fvlFZWZm3tqqqKpqbm3Nr/jSy9+zfs29/Zs2aFRUVFbnb4MGDD3dsAAAAOKoOO7Tr6+vjV7/6VTz88MMp59mv6dOnR2tra+62ZcuWo/49AQAA4HAc1t/RnjJlSixevDieffbZOOWUU3Lbq6uro6OjI3bs2JF3VbulpSWqq6tza55//vm8x9vzqeR71uytrKwsysrKDmdUAAAA6FIFXdHOZrMxZcqUeOyxx+Lpp5+OYcOG5e0fOXJk9OnTJxobG3PbNm7cGJs3b45MJhMREZlMJtatWxfbt2/PrVm2bFmUl5fH8OHDj+S5AAAAQNEVdEW7vr4+Fi1aFD/72c/i+OOPz72nuqKiIo455pioqKiISZMmxbRp02LAgAFRXl4e119/fWQymRg9enRERIwdOzaGDx8eV111VcyePTuam5tjxowZUV9f76o1AAAAPV5BoT1v3ryIiPjMZz6Tt33BggVxzTXXRETEPffcE7169Yrx48dHe3t71NbWxty5c3Nre/fuHYsXL47JkydHJpOJ/v37x8SJE2PmzJlH9kwAAACgGygotLPZ7CHX9OvXLxoaGqKhoeGAa4YOHRpPPPFEId8aAAAAeoQj+jvaAAAAQD6hDQAAAAkJbQAAAEhIaAMAAEBCQhsAAAASEtoAAACQkNAGAACAhIQ2AAAAJCS0AQAAICGhDQAAAAkJbQAAAEhIaAMAAEBCpcUeAOg5Tr1pSbFHyPO7u+qKPQIAAOzDFW0AAABISGgDAABAQkIbAAAAEhLaAAAAkJDQBgAAgISENgAAACQktAEAACAhoQ0AAAAJCW0AAABISGgDAABAQqXFHgDoOaaWPlrsEfZSV+wBAABgH65oAwAAQEJCGwAAABIS2gAAAJCQ0AYAAICEhDYAAAAkJLQBAAAgIaENAAAACQltAAAASEhoAwAAQEJCGwAAABIS2gAAAJCQ0AYAAICEhDYAAAAkJLQBAAAgIaENAAAACQltAAAASEhoAwAAQEJCGwAAABIS2gAAAJCQ0AYAAICEhDYAAAAkJLQBAAAgIaENAAAACQltAAAASEhoAwAAQEJCGwAAABIS2gAAAJCQ0AYAAICEhDYAAAAkJLQBAAAgodJiDwAAAD3BnBlfKvYIeaZ+c0GxRwAOwBVtAAAASEhoAwAAQEIFh/azzz4bF198cQwaNChKSkri8ccfz9ufzWbjlltuiZNPPjmOOeaYqKmpiVdeeSVvzeuvvx4TJkyI8vLyqKysjEmTJsWbb755RE8EAAAAuoOCQ/utt96Ks88+OxoaGva7f/bs2XHvvffG/fffH2vWrIn+/ftHbW1t7Ny5M7dmwoQJsX79+li2bFksXrw4nn322bjuuusO/1kAAABAN1Hwh6GNGzcuxo0bt9992Ww25syZEzNmzIhLLrkkIiJ+8IMfRFVVVTz++ONxxRVXxK9//etYunRpvPDCC3HuuedGRMR9990Xn/3sZ+M73/lODBo06AieDgAAABRX0vdob9q0KZqbm6Ompia3raKiIkaNGhVNTU0REdHU1BSVlZW5yI6IqKmpiV69esWaNWv2+7jt7e3R1taWdwMAAIDuKGloNzc3R0REVVVV3vaqqqrcvubm5hg4cGDe/tLS0hgwYEBuzd5mzZoVFRUVudvgwYNTjg0AAADJ9IhPHZ8+fXq0trbmblu2bCn2SAAAALBfSUO7uro6IiJaWlrytre0tOT2VVdXx/bt2/P2v/322/H666/n1uytrKwsysvL824AAADQHSUN7WHDhkV1dXU0NjbmtrW1tcWaNWsik8lEREQmk4kdO3bE2rVrc2uefvrp6OzsjFGjRqUcBwAAALpcwZ86/uabb8arr76a+3rTpk3x8ssvx4ABA2LIkCExderU+OY3vxmnnXZaDBs2LG6++eYYNGhQXHrppRERceaZZ8ZFF10U1157bdx///2xa9eumDJlSlxxxRU+cRwAAIAer+DQ/uUvfxkXXHBB7utp06ZFRMTEiRNj4cKF8Y1vfCPeeuutuO6662LHjh3xyU9+MpYuXRr9+vXL3eeHP/xhTJkyJS688MLo1atXjB8/Pu69994ETwcAAACKq+DQ/sxnPhPZbPaA+0tKSmLmzJkxc+bMA64ZMGBALFq0qNBvDQAAAN1ej/jUcQAAAOgpCr6iDSSyYlaxJ8h3wfRiTwAAAO8KQhuKZE7jb4o9Qp6pFxx6DQAAcGheOg4AAAAJCW0AAABISGgDAABAQkIbAAAAEhLaAAAAkJDQBgAAgISENgAAACQktAEAACAhoQ0AAAAJCW0AAABISGgDAABAQkIbAAAAEhLaAAAAkJDQBgAAgISENgAAACQktAEAACAhoQ0AAAAJCW0AAABISGgDAABAQkIbAAAAEhLaAAAAkFBpsQcAAACOkhWzij1BvgumF3sC6BKuaAMAAEBCQhsAAAAS8tJxAAB4l5rT+Jtij5Bn6gXFngC6hivaAAAAkJDQBgAAgISENgAAACQktAEAACAhoQ0AAAAJCW0AAABISGgDAABAQv6ONu8Kp960pNgj5PndXXXFHgEAACgSV7QBAAAgIaENAAAACQltAAAASMh7tHlXmFr6aLFH2Iv3aAMAwHuVK9oAAACQkNAGAACAhLx0HAAA6D5WzCr2BPkumF7sCeiBXNEGAACAhIQ2AAAAJCS0AQAAICGhDQAAAAn5MDQAAKDbmNP4m2KPkGfqBcWegJ5IaLMvn/QIAABw2Lx0HAAAABIS2gAAAJCQl46zD++LAQAAOHyuaAMAAEBCQhsAAAASEtoAAACQkPdoAwAAHAl/Hpe9CG0AAIAj4MOE2VtRXzre0NAQp556avTr1y9GjRoVzz//fDHHAQAAgCNWtCvaP/7xj2PatGlx//33x6hRo2LOnDlRW1sbGzdujIEDBxZrrOROvWlJsUfI87u76oo9AgAAwLta0UL77rvvjmuvvTa+9KUvRUTE/fffH0uWLIkHHnggbrrppmKNldzU0keLPcJehDYAAMDRVJTQ7ujoiLVr18b06f//Jv1evXpFTU1NNDU17bO+vb092tvbc1+3trZGRERbW9vRH/YI7WzvKPYIed7JfzMzHzkzdw0zdw0zdw0zdw0zdw0zdw0zd41368xz75jcBZO8c39787xij3BIe/67ZrPZQ64tyb6TVYlt3bo13v/+98eqVasik8nktn/jG9+IlStXxpo1a/LW33bbbXH77bd39ZgAAACQZ8uWLXHKKaccdE2P+NTx6dOnx7Rp03Jfd3Z2xuuvvx4nnHBClJSUFHGyrtHW1haDBw+OLVu2RHl5ebHH4TA4hj2b49fzOYY9m+PXszl+PZ9j2LM5fulks9l44403YtCgQYdcW5TQPvHEE6N3797R0tKSt72lpSWqq6v3WV9WVhZlZWV52yorK4/miN1SeXm5k6OHcwx7Nsev53MMezbHr2dz/Ho+x7Bnc/zSqKioeEfrivLnvfr27RsjR46MxsbG3LbOzs5obGzMeyk5AAAA9DRFe+n4tGnTYuLEiXHuuefGxz/+8ZgzZ0689dZbuU8hBwAAgJ6oaKF9+eWXx3/+53/GLbfcEs3NzXHOOefE0qVLo6qqqlgjdVtlZWVx66237vPyeXoOx7Bnc/x6PsewZ3P8ejbHr+dzDHs2x684ivKp4wAAAPBuVZT3aAMAAMC7ldAGAACAhIQ2AAAAJCS0AQAAICGhDQAAAAkJ7W6ioaEhTj311OjXr1+MGjUqnn/++YOuf+SRR+KMM86Ifv36xYgRI+KJJ57ooknZ26xZs+K8886L448/PgYOHBiXXnppbNy48aD3WbhwYZSUlOTd+vXr10UT86duu+22fY7FGWeccdD7OP+6l1NPPXWfY1hSUhL19fX7Xe/8K65nn302Lr744hg0aFCUlJTE448/nrc/m83GLbfcEieffHIcc8wxUVNTE6+88sohH7fQ36McvoMdw127dsWNN94YI0aMiP79+8egQYPi6quvjq1btx70MQ/nZzGH51Dn4DXXXLPPsbjooosO+bjOwa5zqGO4v9+JJSUl8e1vf/uAj+kcTE9odwM//vGPY9q0aXHrrbfGiy++GGeffXbU1tbG9u3b97t+1apVceWVV8akSZPipZdeiksvvTQuvfTS+NWvftXFkxMRsXLlyqivr4/Vq1fHsmXLYteuXTF27Nh46623Dnq/8vLy2LZtW+722muvddHE7O3DH/5w3rH4l3/5lwOudf51Py+88ELe8Vu2bFlERHzhC1844H2cf8Xz1ltvxdlnnx0NDQ373T979uy499574/777481a9ZE//79o7a2Nnbu3HnAxyz09yhH5mDH8I9//GO8+OKLcfPNN8eLL74YP/3pT2Pjxo3x+c9//pCPW8jPYg7foc7BiIiLLroo71j86Ec/OuhjOge71qGO4Z8eu23btsUDDzwQJSUlMX78+IM+rnMwsSxF9/GPfzxbX1+f+3r37t3ZQYMGZWfNmrXf9V/84hezdXV1edtGjRqV/eu//uujOifvzPbt27MRkV25cuUB1yxYsCBbUVHRdUNxQLfeemv27LPPfsfrnX/d31e+8pXsBz/4wWxnZ+d+9zv/uo+IyD722GO5rzs7O7PV1dXZb3/727ltO3bsyJaVlWV/9KMfHfBxCv09Sjp7H8P9ef7557MRkX3ttdcOuKbQn8Wksb/jN3HixOwll1xS0OM4B4vnnZyDl1xySXbMmDEHXeMcTM8V7SLr6OiItWvXRk1NTW5br169oqamJpqamvZ7n6amprz1ERG1tbUHXE/Xam1tjYiIAQMGHHTdm2++GUOHDo3BgwfHJZdcEuvXr++K8diPV155JQYNGhQf+MAHYsKECbF58+YDrnX+dW8dHR3x0EMPxZe//OUoKSk54DrnX/e0adOmaG5uzjvHKioqYtSoUQc8xw7n9yhdq7W1NUpKSqKysvKg6wr5WczR9cwzz8TAgQPj9NNPj8mTJ8cf/vCHA651DnZvLS0tsWTJkpg0adIh1zoH0xLaRfZf//VfsXv37qiqqsrbXlVVFc3Nzfu9T3Nzc0Hr6TqdnZ0xderUOP/88+MjH/nIAdedfvrp8cADD8TPfvazeOihh6KzszM+8YlPxO9///sunJaIiFGjRsXChQtj6dKlMW/evNi0aVN86lOfijfeeGO/651/3dvjjz8eO3bsiGuuueaAa5x/3dee86iQc+xwfo/SdXbu3Bk33nhjXHnllVFeXn7AdYX+LOboueiii+IHP/hBNDY2xre+9a1YuXJljBs3Lnbv3r3f9c7B7u3BBx+M448/Pi677LKDrnMOplda7AHg3aS+vj5+9atfHfI9LZlMJjKZTO7rT3ziE3HmmWfGd7/73bjjjjuO9pj8iXHjxuX+fdZZZ8WoUaNi6NCh8ZOf/OQd/d9fupf58+fHuHHjYtCgQQdc4/yDrrFr16744he/GNlsNubNm3fQtX4Wdx9XXHFF7t8jRoyIs846Kz74wQ/GM888ExdeeGERJ+NwPPDAAzFhwoRDfuinczA9V7SL7MQTT4zevXtHS0tL3vaWlpaorq7e732qq6sLWk/XmDJlSixevDhWrFgRp5xySkH37dOnT3z0ox+NV1999ShNxztVWVkZH/rQhw54LJx/3ddrr70Wy5cvj7/6q78q6H7Ov+5jz3lUyDl2OL9HOfr2RPZrr70Wy5YtO+jV7P051M9ius4HPvCBOPHEEw94LJyD3dc///M/x8aNGwv+vRjhHExBaBdZ3759Y+TIkdHY2Jjb1tnZGY2NjXlXXP5UJpPJWx8RsWzZsgOu5+jKZrMxZcqUeOyxx+Lpp5+OYcOGFfwYu3fvjnXr1sXJJ598FCakEG+++Wb89re/PeCxcP51XwsWLIiBAwdGXV1dQfdz/nUfw4YNi+rq6rxzrK2tLdasWXPAc+xwfo9ydO2J7FdeeSWWL18eJ5xwQsGPcaifxXSd3//+9/GHP/zhgMfCOdh9zZ8/P0aOHBlnn312wfd1DiZQ7E9jI5t9+OGHs2VlZdmFCxdmN2zYkL3uuuuylZWV2ebm5mw2m81eddVV2Ztuuim3/rnnnsuWlpZmv/Od72R//etfZ2+99dZsnz59suvWrSvWU3hPmzx5craioiL7zDPPZLdt25a7/fGPf8yt2fsY3n777dknn3wy+9vf/ja7du3a7BVXXJHt169fdv369cV4Cu9pX/3qV7PPPPNMdtOmTdnnnnsuW1NTkz3xxBOz27dvz2azzr+eYvfu3dkhQ4Zkb7zxxn32Of+6lzfeeCP70ksvZV966aVsRGTvvvvu7EsvvZT7ROq77rorW1lZmf3Zz36W/dd//dfsJZdckh02bFj2f/7nf3KPMWbMmOx9992X+/pQv0dJ62DHsKOjI/v5z38+e8opp2RffvnlvN+L7e3tucfY+xge6mcx6Rzs+L3xxhvZr33ta9mmpqbspk2bssuXL89+7GMfy5522mnZnTt35h7DOVhch/o5ms1ms62trdljjz02O2/evP0+hnPw6BPa3cR9992XHTJkSLZv377Zj3/849nVq1fn9n3605/OTpw4MW/9T37yk+yHPvShbN++fbMf/vCHs0uWLOniidkjIvZ7W7BgQW7N3sdw6tSpueNdVVWV/exnP5t98cUXu354spdffnn25JNPzvbt2zf7/ve/P3v55ZdnX3311dx+51/P8OSTT2YjIrtx48Z99jn/upcVK1bs92fmnmPU2dmZvfnmm7NVVVXZsrKy7IUXXrjPcR06dGj21ltvzdt2sN+jpHWwY7hp06YD/l5csWJF7jH2PoaH+llMOgc7fn/84x+zY8eOzZ500knZPn36ZIcOHZq99tpr9wlm52BxHernaDabzX73u9/NHnPMMdkdO3bs9zGcg0dfSTabzR7VS+YAAADwHuI92gAAAJCQ0AYAAICEhDYAAAAkJLQBAAAgIaENAAAACQltAAAASEhoAwAAQEJCGwAAABIS2gAAAJCQ0AYAAICEhDYAAAAk9L8tyaK9Lr5FNQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 1200x700 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize = (12, 7))\n",
    "ax.bar(range(len(expected)), expected)\n",
    "ax.bar(range(len(observed)), observed, alpha = .5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Гипотеза о нормальности выборки уже не отвергается. (Заметьте, мы знаем точно, что она не нормальная). Вот у нас и произошла ошибка первого рода). Тем не менее, на такой выборке мы весьма вероятно сможем применить t-test (который, на самом деле, к некоторым отклонениям от нормальности не так чувствителен)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Проблема вышеописанного метода Хи-квадрат заключается в том, что количество бакетов заранее не определено и при его изменении мы точно будем получать разные результаты. Кроме того, выборка должна быть достаточно большой (формально - в 80% бакетов ожидаемое кол-во наблюдений должно быть не менее чем 5). Но есть более точный подход\n",
    "http://medstatistic.ru/articles/doveritelnye-intervaly-dlya-chastot-i-doley.pdf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## QQ-plot и Критерий Шапиро-Уилка"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Есть намного более удобный способ визуальной и численной проверки нормальности. Это так называемый QQ-график."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "smpl = np.array([pylab.normal(70, 5) + pylab.exponential(1) for _ in range(100)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjIAAAHHCAYAAACle7JuAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABfPklEQVR4nO3de3zO9f/H8ce1sTluDpsZm/MhpZyTUyiFkvNZmZxKzql+JMdCIoRyKoYcY+ikRNRyJodvhRzmNHMIM8eN7fP749OuzDaua67t2rU977fbdeNzuD7Xy9X3a0/vo8UwDAMRERERF+Tm7AJEREREUkpBRkRERFyWgoyIiIi4LAUZERERcVkKMiIiIuKyFGRERETEZSnIiIiIiMtSkBERERGXpSAjIiIiLktBRkRSxGKx0KdPH4c9Lzg4GIvFwq5dux54b7169ahXr571+Pjx41gsFoKDg63nRo4cicVicVh9DyupGkXk4SnIiGQg8WEg/pUtWzbKlClDnz59OHfunLPLc7qxY8eyevVqhz5z06ZNCb7zrFmzUqJECTp37syxY8cc8hlbtmxh5MiRREZGOuR5IhmJgoxIBjR69GgWLlzI9OnTqVmzJjNmzKBGjRrcuHHD2aU5xLp161i3bt1973nvvfe4efNmgnOpEWTi9evXj4ULFzJ79mxefPFFli1bRrVq1Thz5sxDP3vLli2MGjVKQUYkCVmcXYCIOF7jxo2pWrUqAN27dyd//vxMmjSJNWvW0KFDhyTfc/36dXLmzJmWZaaYh4fHA+/JkiULWbKk3V9xderUoXXr1gC8+uqrlClThn79+jF//nyGDBmSZnWIZDZqkRHJBJ555hkAwsLCAOjSpQu5cuXi6NGjvPDCC+TOnZtOnToBZqAZNGgQgYGBeHp6UrZsWSZOnIhhGEk+e9GiRZQtW5Zs2bJRpUoVfv311wTXT5w4wRtvvEHZsmXJnj07+fPnp02bNhw/fjzJ5924cYPXXnuN/Pnz4+XlRefOnbl8+XKCe+4dI5OUe8fIWCwWrl+/zvz5863dQF26dGHjxo1YLBZWrVqV6BmLFy/GYrGwdevW+35WUu79zpPz888/U6dOHXLmzEmePHlo1qwZBw4cSPDnePvttwEoXry4tfbkvj+RzEYtMiKZwNGjRwHInz+/9dydO3do2LAhtWvXZuLEieTIkQPDMGjatCkbN26kW7duVKxYkR9//JG3336b8PBwJk+enOC5v/zyC8uWLaNfv354enry2Wef0ahRI3bs2EH58uUB2LlzJ1u2bKF9+/YEBARw/PhxZsyYQb169fjrr7/IkSNHgmf26dOHPHnyMHLkSA4dOsSMGTM4ceKEdSxKSi1cuJDu3bvz5JNP0rNnTwBKlizJU089RWBgIIsWLaJFixYJ3rNo0SJKlixJjRo17P68pL7ze61fv57GjRtTokQJRo4cyc2bN5k2bRq1atXi999/p1ixYrRs2ZK///6bJUuWMHnyZHx8fADw9fW1uyaRDMkQkQxj3rx5BmCsX7/euHDhgnHq1Clj6dKlRv78+Y3s2bMbp0+fNgzDMIKCggzAGDx4cIL3r1692gCMDz74IMH51q1bGxaLxThy5Ij1HGAAxq5du6znTpw4YWTLls1o0aKF9dyNGzcS1bl161YDMBYsWJCo9ipVqhgxMTHW8x999JEBGGvWrLGeq1u3rlG3bl3rcVhYmAEY8+bNs54bMWKEce9fcTlz5jSCgoIS1TNkyBDD09PTiIyMtJ47f/68kSVLFmPEiBGJ7r/bxo0bDcCYO3euceHCBePMmTPGd999ZxQrVsywWCzGzp07k62xYsWKRoECBYyLFy9az+3bt89wc3MzOnfubD03YcIEAzDCwsLuW4tIZqSuJZEMqEGDBvj6+hIYGEj79u3JlSsXq1atonDhwgnu69WrV4Lj77//Hnd3d/r165fg/KBBgzAMg7Vr1yY4X6NGDapUqWI9LlKkCM2aNePHH38kNjYWgOzZs1uv3759m4sXL1KqVCny5MnD77//nqj2nj17kjVr1gQ1ZsmShe+//97Ob8F2nTt3Jjo6mhUrVljPLVu2jDt37vDyyy/b9IyuXbvi6+tLoUKFePHFF63dWPFjle4VERHB3r176dKlC/ny5bOef+KJJ3juuedS9c8rkpGoa0kkA/r0008pU6YMWbJkwc/Pj7Jly+LmlvDfLVmyZCEgICDBuRMnTlCoUCFy586d4Hy5cuWs1+9WunTpRJ9dpkwZbty4wYULFyhYsCA3b95k3LhxzJs3j/Dw8ARjba5cuZLo/fc+M1euXPj7+6fqmJBHHnmEatWqsWjRIrp16waY3UpPPfUUpUqVsukZw4cPp06dOri7u+Pj40O5cuXuO9g4/rssW7ZsomvlypXjxx9/dKkB2CLOoiAjkgE9+eSTybYExPP09EwUblJD3759mTdvHgMGDKBGjRp4e3tjsVho3749cXFxqf75turcuTP9+/fn9OnTREdHs23bNqZPn27z+x9//HEaNGiQihWKSFIUZETEqmjRoqxfv56rV68maJU5ePCg9frdDh8+nOgZf//9Nzly5LAORl2xYgVBQUF8/PHH1ntu3bqV7Joohw8fpn79+tbja9euERERwQsvvJDiP1e8+w0Wbt++PW+++SZLlizh5s2bZM2alXbt2j30ZyYn/rs8dOhQomsHDx7Ex8fH2hqTnlYoFklvNEZGRKxeeOEFYmNjE7VETJ48GYvFQuPGjROc37p1a4JxLqdOnWLNmjU8//zzuLu7A+Du7p5o6va0adOsY2juNXv2bG7fvm09njFjBnfu3En02SmRM2fOZAOUj48PjRs35ssvv2TRokU0atTIOkMoNfj7+1OxYkXmz5+foKY//viDdevWJQhu8YFGC+KJJKYWGRGxeumll6hfvz5Dhw7l+PHjVKhQgXXr1rFmzRoGDBhAyZIlE9xfvnx5GjZsmGD6NcCoUaOs9zRp0oSFCxfi7e3No48+ytatW1m/fn2y05JjYmJ49tlnadu2LYcOHeKzzz6jdu3aNG3a9KH/fFWqVGH9+vVMmjSJQoUKUbx4capXr2693rlzZ+uidu+///5Df96DTJgwgcaNG1OjRg26detmnX7t7e3NyJEjE9QNMHToUNq3b0/WrFl56aWXNH5GBAUZEbmLm5sbX3/9NcOHD2fZsmXMmzePYsWKMWHCBAYNGpTo/rp161KjRg1GjRrFyZMnefTRRwkODuaJJ56w3vPJJ5/g7u7OokWLuHXrFrVq1WL9+vU0bNgwyRqmT5/OokWLGD58OLdv36ZDhw5MnTrVId0rkyZNomfPntbtC4KCghIEmZdeeom8efMSFxfnkOD0IA0aNOCHH35gxIgRDB8+nKxZs1K3bl3Gjx9P8eLFrfdVq1aN999/n5kzZ/LDDz8QFxdHWFiYgowIYDHubfMVEcmk7ty5Q6FChXjppZf44osvnF2OiNhAY2RERP61evVqLly4QOfOnZ1diojYSC0yIpLpbd++nf379/P+++/j4+OT5EJ9IpI+qUVGRDK9GTNm0KtXLwoUKMCCBQucXY6I2EEtMiIiIuKy1CIjIiIiLktBRkRERFxWhl9HJi4ujjNnzpA7d24t8y0iIuIiDMPg6tWrFCpU6L77wmX4IHPmzBkCAwOdXYaIiIikwKlTpwgICEj2eoYPMvEb3506dQovLy8nVyMiIiK2iIqKIjAwMMEGtknJ8EEmvjvJy8tLQUZERMTFPGhYiAb7ioiIiMtSkBERERGXpSAjIiIiLktBRkRERFyWgoyIiIi4LAUZERERcVkKMiIiIuKyFGRERETEZSnIiIiIiMtSkBERERGXpSAjIiIiLktBRkRERFyWgoyIiIikTGwsfPedU0tQkBERERH7nT0LDRtCkyawdKnTysjitE8WERER17RhA3TqBOfOQY4cEBfntFLUIiMiIiK2iY2FESPguefMEFO+POzcCR07Oq0ktciIiIjIg505YwaWX34xj7t3h08+MVtknEhBRkRERO7vxx/h5Zfhn38gVy6YNcuprTB3U5ARERGRpN25A8OGwYcfmscVKsDy5VCmDLGxEBoKERHg7w916oC7e9qXqCAjIiIiiZ06BR06wObN5nGvXjBpEmTLRkgI9O8Pp0//d3tAgNnT1LJl2papwb4iIiKS0HffQcWKZojx8oJly+Czz6whpnXrhCEGIDzcPB8SkralKsiIiIiI6fZteOstc22YS5egShX4/Xdo2xYwJy317w+Gkfit8ecGDDDvSysKMiIiIgLHj5sDXT7+2Dzu189skSlZ0npLaGjilpi7GYbZIxUamrql3k1jZERERDK71avh1VchMhLy5IF586B580S3RUTY9jhb73MEtciIiIhkVtHRZl9RixZmiKleHfbsSTLEgDk7yRa23ucICjIiIiKZ0dGjUKsWTJ1qHg8aBL/+CsWKJfuWOnXM2UkWS9LXLRYIDDTvSysKMiIiIpnNV19B5cqwezfkywfffAMTJ4KHx33f5u5uTrGGxGEm/njKlLRdT0ZBRkREJLO4dQveeMOchRQVZbbI7N1rzlKyUcuWsGIFFC6c8HxAgHk+rdeR0WBfERGRzODwYTPA7N1rHg8ZAqNGQdasdj+qZUto1kwr+4qIiEhaWLIEevaEa9fAxwe+/BIaNnyoR7q7Q716jinvYahrSUREJKO6cQN69DA3eLx2DerWhX37HjrEpCcKMiIiIhnRgQPmdOrPPzdH4g4fDuvXQ6FCzq7ModS1JCIiktHMn28O6r1xA/z8YNEiePZZZ1eVKtQiIyIiklFcvw5dupivGzfM8LJ3b4YNMaAgIyIikjH88QdUq2a2xri5wejR8OOPULCgsytLVepaEhERcWWGAXPnQp8+5joxhQrB4sXmwN5MQEFGRETEVV29Cr16mWNgABo1ggULwNfXuXWlIXUtiYiIuKJ9+6BqVTPEuLvDuHHw3XeZKsSAWmRERERci2HArFkwYIC5e3VAACxdam43kAkpyIiIiLiKK1fMFXqXLzePmzSB4GDIn9+pZTmTupZERERcwe7dUKWKGWKyZIGPP4avv87UIQbUIiMiIpK+GQZMnw5vvQUxMVC0KCxbZq7aKwoyIiIi6dbly9CtG6xaZR43b25Otc6b16llpSfqWhIREUmPduyAypXNEOPhAVOnQkiIQsw9FGRERETSE8OASZPMWUjHj0OJErBlC/Tta27+KAmoa0lERCS9uHTJ3Cfpm2/M4zZtYM4c8PZ2alnpmVpkRERE0oMtW6BiRTPEeHrCjBnmoF6FmPtSkBEREXGmuDgYPx6efhpOnYLSpWHbNnj9dXUl2UBdSyIiIs5y4QIEBcHateZxx44wcybkzu3culyIgoyIiIgz/PordOgAZ85AtmzmWjFdu9rdChMbC6GhEBEB/v5Qp4659VJmoa4lERGRtBQbCx98APXrmyHmkUdg505zvRg7Q0xICBQrZj6qY0fz12LFzPOZhYKMiIhIWjl3Dho1gmHDzLExQUGwaxeUL2/3o0JCoHVrOH064fnwcPN8ZgkzCjIiIiJp4eefzVlJ69dDjhzmZo/BwZAzp92Pio2F/v3NJWfuFX9uwADzvoxOQUZERCQ1xcbCiBHQoAGcPWu2vuzcabbGpFBoaOKWmLsZhjkBKjQ0xR/hMjTYV0REJLWcOQOdOsGmTeZx9+7wySdmi8xDiIhw7H2uTEFGREQkNaxbBy+/bE6xzpULZs0yR+Q6gL+/Y+9zZepaEhERcaQ7d+Ddd6FhQzPEVKgAu3c7LMSAOcU6ICD5SU4WCwQGmvdldE4NMrGxsQwbNozixYuTPXt2SpYsyfvvv49x1+glwzAYPnw4/v7+ZM+enQYNGnD48GEnVi0iIpKM06fNOdDjxpnHr79urtJbpkyiW2NjzR6nJUvMX+0ZmOvubvZQQeIwE388ZUrmWE/GqUFm/PjxzJgxg+nTp3PgwAHGjx/PRx99xLRp06z3fPTRR0ydOpWZM2eyfft2cubMScOGDbl165YTKxcREbnHd9+Zs5J++81cmXfZMnO/pGzZEt3qiPVfWraEFSugcOGE5wMCzPMtWz7MH8Z1WAwjqclbaaNJkyb4+fnxxRdfWM+1atWK7Nmz8+WXX2IYBoUKFWLQoEG89dZbAFy5cgU/Pz+Cg4Np3779Az8jKioKb29vrly5gpeXV6r9WUREJJO6fdvsSpo40TyuXBmWL4eSJZO8PX79l3t/+sa3pNgbQjLqyr62/vx2aotMzZo12bBhA3///TcA+/bt47fffqNx48YAhIWFcfbsWRo0aGB9j7e3N9WrV2fr1q1OqVlERMTqxAkzOcSHmL59zV2skwkxqbH+i7s71Ktn7nZQr17GCDH2cOqspcGDBxMVFcUjjzyCu7s7sbGxjBkzhk6dOgFw9uxZAPz8/BK8z8/Pz3rtXtHR0URHR1uPo6KiUql6ERHJ1FavhldfhchIyJMH5s6FFi3u+xZ71n+pV8+BtWZgTm2RWb58OYsWLWLx4sX8/vvvzJ8/n4kTJzJ//vwUP3PcuHF4e3tbX4GBgQ6sWEREMr2YGLPZpEULM8Q8+STs2fPAEANa/yU1ODXIvP322wwePJj27dvz+OOP88orrzBw4EDG/Tvau2DBggCcO3cuwfvOnTtnvXavIUOGcOXKFevr1KlTqfuHEBGRzOPYMahV678pQ4MGmc0nxYrZ9Hat/+J4Tg0yN27cwM0tYQnu7u7ExcUBULx4cQoWLMiGDRus16Oioti+fTs1atRI8pmenp54eXkleImIiDy0FSugUiVzk8d8+eDrr82xMR4eNj9C6784nlPHyLz00kuMGTOGIkWK8Nhjj7Fnzx4mTZpE165dAbBYLAwYMIAPPviA0qVLU7x4cYYNG0ahQoVo3ry5M0sXEZHM4tYts+Xls8/M45o1YelSM3HYKX79l9atzdBy96DfzLb+i6M4dfr11atXGTZsGKtWreL8+fMUKlSIDh06MHz4cDz+TbiGYTBixAhmz55NZGQktWvX5rPPPqNMEosLJUXTr0VEJMUOH4a2bWHvXvN48GAYPRqyZk1wm71ToENCzNlLdw/8DQw0Q0xmWf/lQWz9+e3UIJMWFGRERCRFliyBnj3h2jXw8YGFC6FRo0S3JRVKAgLMlpf7hZKMuv6LoyjI/EtBRkRE7HLzpplM5swxj+vWhcWLoVChRLc6enE7+Y9LLIgnIiKSrhw8aE6nnjPHTCPDhsH69UmGmNRY3E7spyAjIiICsGABVKkCf/wBfn6wbp05HiZL0vNi7FncTlKPgoyIiGRu16+bK/QGBcGNG/DMM+bg3ru2x0mKFrdLHxRkREQk8/rzT6hWDYKDwc3NbIFZtw6SWXT1blrcLn1w6joyIiIiTmEY5t5Iffuag3v9/c0BvfdscHS/mUXxi9uFhyc9TsZiMa9rcbvUpRYZERHJXK5ehVdege7dzRDz/PNmV9I9ISYkxNx5oH596NjR/LVYMfM8/Le4HSReqVeL26UdBRkREck89u2DqlVh0SIzYYwdC2vXQoECCW6Ln1Z972De8HDzfHyYadnSnGJduHDC+wICNPU6rWgdGRERyfgMA2bPNudLR0ebyWPpUqhdO9GtsbFmy0tyM5Liu4zCwv5rbdHido5n689vjZEREZGMLSoKevSA5cvN4xdfNAf3+vgkebs906rje6Pc3RP1TEkaUdeSiIhkXL//DpUrmyEmSxaYMMHctTqZEAOaVu1q1CIjIiIZj2HA9Onw1lsQEwNFi5pdSU899cC3alq1a1GQERGRjCUyErp1+29EbrNmMG8e5M2b4LbkxrVoWrVrUdeSiIhkHDt2QKVKZojJmtWc/7xqVaIQc7+p1ZpW7VoUZERExPUZBkyaBLVqwfHjULw4sb9uZlOF/ixZamHTpv82b7RlarWmVbsOTb8WERHXdukSdOkC33xjHrduzddNP6f3u94JwkpAgJl13nzT9qnVmlbtPLb+/FaQERER17VlC7Rvb86H9vCAyZMJ8etF6zaWRONbLJakx7wkZeNGTad2Nq0jIyIiGVdcnDmVeuhQs9mkVClYvpzYJyrRv1jSgcWef7ZrarXrUJARERHXcuECBAWZWwsAce07sPmVWZw+mJtzv9x/MTtbaWq161CQERER1/Hrr9ChA5w5A9my8XuXqTT7pjunl1oe/F4baGq161GQERGRdCXJAbbEEjdmHJZRI7DExXEjsCyb+y+n4dtP2NVldLd7x8xoarVrUpAREZE0kVRAgYTn/vkHBg5M2D1U0f8cyzxepsyJ9QAs4BXeOPUZt/4vl90hJr7FZdKkxJ8TEGCGGE2tdi0KMiIikiruDi6HD8OcOQmDQ/785q8XLyb/jPr8zOKIjhTkHDfITm8+JZhX//0A++q5u8WlZUto0UJTqzMCBRkREXG4kBDo3//+A2/vF2DciGU4oxnG+7hh8AeP0ZblHODRFNd0b4uLdqzOGBRkRETEoeJXzk3p2BV/zrCITtRnEwBf0JW+TOMmOex+1uTJ4OenFpeMTEFGREQcJjbWbIlJaYh5jnV8ycsU4ALXyMnrzGQRL9v9nPixMH37KrxkdNprSUREHCY0NGXruLhzhw8Yyg80ogAX2McTVGF3ikMMaPZRZqEgIyIiDpOSFXELc5qN1GcoY3HDYCav8RTb+JuyNr3/3rCijR0zF3UtiYiI3e6dSl2zprnt0V9/2fecxnzPAjrjw0WiyE0P5rCcdja9N77lZckS8PXV7KPMSkFGRETsktSMpPidom2VhduMYSjvMAGA3VSmHcs4SqlE9wYGmvtCLlmidV8kMQUZERGxWXIzkuwJMUU4wVLaU4NtAEylL28zgRg8CQyEjz9OuoVl3Dit+yKJKciIiIhNHnZGUkAATKizhhZfv4rn9csY3t7EfT6XJ3xaEmxDONG6L5IUBRkREbFJSmckvfceNHg6hjrfvoPb1E/Mk9WqYVm2DPfixann0Cols1GQERGRJN09oLdAAfj555Q9p5pPGHWHtoOdO80TAwfChx+Ch4fjipVMS0FGREQSsWWLAVu0ZCUvDO0G169A3rwQHAxNmzqkRhFQkBERkX/Ft8CsWWPOBnoYntziY96iN5/CdaBGDVi6FIoUcUSpIlYKMiIi4rAWGIBSHGYZ7ajMHvPEO+/ABx9A1qwP/3CReyjIiIhkcg+7ySP8t45MO5Yyhx7k5hrRXj54Ll0AjRs7rliRe2iLAhGRTOxhp1SDOSvpxsWbhDd5jaV0IDfXMGrXwfOvvQoxkuoUZEREMrGUTqm+20ulD+JRpzqFvp1t7hswdCiWjT9D4cKOKVLkPtS1JCKSiYWHp/y9Fgv0y7uQam/0guvXzTnaCxfC8887rkCRB1CQERHJpEJCYMCAlL03BzeYZvSl66W55on69WHRInN5XpE0pCAjIpIJPcwA30f5k1VZ2lLmzl/g5gbDh5sDZbTxkTiBgoyISCZj7wDfAQOgSRPAMMj5VTDV5vfGPfomFCxobkmtDZDEiRRkREQyGVsH+Pr6wsyZ0LIlcO0a9OoFX35pXnz+eXM8TIECqVqryINo1pKISCYSGwsbNth27+TJ/4aY/fuhShUzxLi5wZgxsHatQoykC2qRERHJJOxdvbdwIQNmzzHfdOuWOZ16yRKoUyd1CxWxg4KMiEgmYM/gXosFHikUxdOzXoNlS82TL7wA8+eDj0/qFipiJ3UtiYhkcPYM7rVYoJLxOztiq+C2bClkyQITJsA33yjESLqkICMiksHZvnqvwbven7IjSw1ynT1i7lQdGgpvvWWOjRFJh9S1JCKSwUVEPPgebyLZ/Eh3Hju40jzRrBnMnQv58qVucSIPSRFbRCSDe9Biu1XZye9UNkNM1qzmdKVVqxRixCUoyIiIZHB16kBAgDn+JSGD/kxhM7UoQRhG8eKwebO5Al7im0XSJQUZEZEMzt0dPvnE/H18PsnLJVbTnCkMxIPbhD/VCsvvv0O1as4rVCQFFGRERDK42Fizl6h/f3Pi0VNsZQ+VaMbXROPB3u7TKbzlK8iTx9mlithNg31FRDKg2FhzwtGaNeam1BcugIU43mIiY3mXLMRyo3ApPFcvp2LVSs4uVyTFFGRERDKYpFbwzc8/LKAzL7AWgKW0I8eHs2la1ctJVYo4hrqWREQyiNhYGD0aWrVKGGJqE8peKvICa7lJNnoyi44soc+7XsTGOq9eEUdwapApVqwYFosl0at3794A1KtXL9G1119/3Zkli4ikK7GxsGkTDBwIBQvCiBH/XbMQx7uMYRP1CCCcg5SlOtuZQ08MLJw6ZXY/ibgyp3Yt7dy5k9i7/jnwxx9/8Nxzz9GmTRvruR49ejB69GjrcY4cOdK0RhGR9Cap8S/3KsA5FvIKz/MTAAt5mV7M4Dq5Etxny2J5IumZU4OMr69vguMPP/yQkiVLUrduXeu5HDlyULBgwbQuTUQkXbJlB+t6bGQxHfHnLDfITm8+JZguQOK1YR60WJ5IepduxsjExMTw5Zdf0rVrVyx3LcS0aNEifHx8KF++PEOGDOHGjRtOrFJExHnid7BOLsS4EctwRrGeBvhzlj95lGrsJJhXuTfEWCwQGGgulifiyuxukTl16hQWi4WAgAAAduzYweLFi3n00Ufp2bNnigtZvXo1kZGRdOnSxXquY8eOFC1alEKFCrF//37+7//+j0OHDhESEpLsc6Kjo4mOjrYeR0VFpbgmEZH0IiYGXn89+R2sCxLBIjrxDBsB+IKu9GUaN0ncHR//b8UpU8zF8kRcmmGn2rVrGwsWLDAMwzAiIiIMLy8vo0aNGoaPj48xatQoex9n9fzzzxtNmjS57z0bNmwwAOPIkSPJ3jNixAgDSPS6cuVKimsTEXGGO3cMY+NGwxgwwDC8vQ3DjDGJXw1YZ5ylgGGAcZWcxsssSPZeMIzAQMNYudLZfzqR+7ty5YpNP78thpFcvk9a3rx52bZtG2XLlmXq1KksW7aMzZs3s27dOl5//XWOHTtmd5g6ceIEJUqUICQkhGbNmiV73/Xr18mVKxc//PADDRs2TPKepFpkAgMDuXLlCl5eWi9BRNKn+AG84eHm4N3jx2Hx4qQH8sZz5w4jGcm7jMUNg308QTuWcYhHEt3r4wMvv2xual2njlpiJP2LiorC29v7gT+/7e5aun37Np6engCsX7+epk2bAvDII48QkcLh7/PmzaNAgQK8+OKL971v7969APjfZ3Sap6entT4REVdgywDeexUinCV04GnM+dOz6MkApnCL7InuHTUKhg5VeJGMye4g89hjjzFz5kxefPFFfvrpJ95//30Azpw5Q/78+e0uIC4ujnnz5hEUFESWLP+Vc/ToURYvXswLL7xA/vz52b9/PwMHDuTpp5/miSeesPtzRETSo/gBvPa0jTdiLQvojC//EEVuejKbZbRPdF9goDkOpmVLx9Urkt7YHWTGjx9PixYtmDBhAkFBQVSoUAGAr7/+mieffNLuAtavX8/Jkyfp2rVrgvMeHh6sX7+eKVOmcP36dQIDA2nVqhXvvfee3Z8hIpIePWgA772ycJsPeI//4yMAdlOZdizjKKWs9/j6QqdO6kKSzMPuMTIAsbGxREVFkTdvXuu548ePkyNHDgoUKODQAh+WrX1sIiJp4e7F7ObNgytXbHtfICdZQgdqsQWAafThLSYSgyfe3vDqqwovkrGk2hgZAMMw2L17N0ePHqVjx47kzp0bDw8PrborInIfKRkLA9CEb5hPEPm4TCTedGUuqzD7i3x9zed5eKRCwSIuwO4gc+LECRo1asTJkyeJjo7mueeeI3fu3IwfP57o6GhmzpyZGnWKiLi0FSvgrt1XbJKVGMYxhEFMAmAH1WjHMo5T3LoWzMyZCjGSudm9sm///v2pWrUqly9fJnv2/0bHt2jRgg0bNji0OBGRjOCrr6B94rG491WMMEKpYw0xkxhIbX7jOMUBCAgww5EG8kpmZ3eLTGhoKFu2bMHjnn8CFCtWjPDwcIcVJiKSEYSEQNu29r2nBSHMpSt5uMIl8tKFYL6hqQbyiiTB7iATFxeXYMfqeKdPnyZ37twOKUpEJCOIn5VkK09uMYG36ct0ALZQg975llKvcxE2KryIJMnurqXnn3+eKVOmWI8tFgvXrl1jxIgRvPDCC46sTUTEZYWEQOHC91+Z924lOcIWalpDzB8vvsPtn35h1/kiTJ4M9eopxIgkxe7p16dPn6Zhw4YYhsHhw4epWrUqhw8fxsfHh19//VXTr0Uk07N3YG87ljKbnnhxlX/Iz6F3F1BrjP5hKJmbrT+/U7SOzJ07d1i6dCn79+/n2rVrVK5cmU6dOiUY/JteKMiISFr66ivo0MFcL+ZBsnGTyQzkdWYBsN2zDpemLaZxj4BUrlIk/UvVdWSyZMnCyy+/nOLiREQyirs3e9ywwVzkzhZlOMRy2lKB/cRh4dTL71L185G4e6bor2WRTMvu/8csWLDgvtc7d+6c4mJERFxJShe468SXzOR1cnGdW16+ZPvqS4o+/3zqFCmSwdndtXT3tgRg7oZ948YN68q+ly5dcmiBD0tdSyKSGlKy2WMOrjONvnTFbLY5X74+BdYtAn//VKpSxHXZ+vPb7llLly9fTvC6du0ahw4donbt2ixZsuShihYRcQWxsWZLjD0hphx/sYMn6co84rAQ+94ICuz9SSFG5CE5pDO2dOnSfPjhh7z88sscPHjQEY8UEUk37h4Hc+ECnD9vT3eSQReC+ZTe5OAmERTkyMhF1BnxTGqWLJJpOGxUWZYsWThz5oyjHiciki6kdBwMQE6u8Rlv0JmFAKynATfnfMlL3f0cXKVI5mV3kPn6668THBuGQUREBNOnT6dWrVoOK0xExNlSMg4m3uPsZzlteYRDxOLGMN6n0tLBtGlnd4++iNyH3UGmefPmCY4tFgu+vr4888wzfPzxx46qS0TEqVIyDsZk0IM5fEJ/snOL0xSmv+8SOs2sow0eRVJBivZaEhHJ6DZtsr87KTdRzOI1OrAUgO9pzJ/vLGD5WB9tLyCSStTGKSJyj5TsWF2RPeymCh1Yyh3cGeP9Ebe++pa3xyvEiKQmm1pk3nzzTZsfOGnSpBQXIyLibPbukwQGvZjBZAbiSQyRXoEcG7uMwa/XUIARSQM2BZk9e/bY9DCLxfJQxYiIOFP8Pkm28iaSOfSgDSsAWJetKc8emUdl33ypVKGI3MumILNx48bUrkNExKns7U6qyk6W0Y4ShBFDVgYzntpfDsDdV/+gE0lLGiMjIple/Awl2xj0ZwqbqUUJwgijGC0LbKb2yoG0bKUQI5LWUrQg3q5du1i+fDknT54kJiYmwbWQkBCHFCYiklZCQ22boZSXS8ylK81ZA8DJai0JH/4Faxrn0XgYESexu0Vm6dKl1KxZkwMHDrBq1Spu377Nn3/+yc8//4y3t3dq1CgikqrCwx98T3W2sYdKNGcNhocHTJ9Oke0rqN1EIUbEmewOMmPHjmXy5Ml88803eHh48Mknn3Dw4EHatm1LkSJFUqNGEZFUERsLo0dDr17J32MhjreYQCh1KMpJrhUsiWXrVujdGzTBQcTp7A4yR48e5cUXXwTAw8OD69evY7FYGDhwILNnz3Z4gSIiqSEkBPz8YMQIuHo16XvycZGvacoE3iErdzhVsx25Dv0OlSunbbEikiy7g0zevHm5+u//6wsXLswff/wBQGRkJDdu3HBsdSIiqSAkBFq1gosXk7+nFr+xl4o04Ttu4cnuHjMJ/G0JeHmlXaEi8kB2B5mnn36an376CYA2bdrQv39/evToQYcOHXj22WcdXqCIiKPExsKGDdCjR/L3WIhjMOPYRD0COc0R9zL89vEOqsx+TV1JIumQzbOW/vjjD8qXL8/06dO5desWAEOHDiVr1qxs2bKFVq1a8d5776VaoSIiDyMkxJxifb/ZSb6cZyGv0JB1ACzgFTxnfUa7brnSqEoRsZfNQeaJJ56gWrVqdO/enfbt2wPg5ubG4MGDU604ERFHCAmB1q3vv5N1XTaxmI4UIoIbZKc3nxJMFzaWVCuMSHpmc9fSL7/8wmOPPcagQYPw9/cnKCiI0NDQ1KxNROShxcTA668nH2LciGUYo9nAsxQigj95lGrsJJhXCQy0UKdO2tYrIvaxOcjUqVOHuXPnEhERwbRp0zh+/Dh169alTJkyjB8/nrNnz6ZmnSIiNomNhU2bYNEi6NoVfHzgwoWk7/XjLOt4ntGMwJ045vIqT7KDv3gMiwWmTEFrxIikcxbDuF9j6/0dOXKEefPmsXDhQs6ePUujRo34+uuvHVnfQ4uKisLb25srV67gpdkGIhmaLeNg4j3LehbRCT/Oc42c9GIGX/IKAPnzw+zZ0LJlKhcsIsmy9ef3QwUZgOvXr7No0SKGDBlCZGQksbGxD/M4h1OQEcm4YmPN7QXCw83ZSPPmPfg97txhJCN5l7G4YbCfx2nLcg7xCABdusDnn6slRsTZbP35naK9lgB+/fVX5s6dy8qVK3Fzc6Nt27Z069YtpY8TEbGLPa0v8QoRzhI68DTm+L5Z9GQAU7hFdgACAhRiRFyNXUHmzJkzBAcHExwczJEjR6hZsyZTp06lbdu25MyZM7VqFBFJYMUKaNPGvvc0Yi0L6Iwv/3CVXPRgDstob71uscAnnyjEiLgam4NM48aNWb9+PT4+PnTu3JmuXbtStmzZ1KxNRCSRr76CDh1svz8Lt3mfYQxmPAB7qEhblnOE0tZ7AgPNgb0aEyPiemwOMlmzZmXFihU0adIEd/2TRUTSWGwsjBlj7o1kq0BOsoQO1GILANPpzVtMJJpseHvDtGlmiKlTRy0xIq7K5iCT3mYjiUjmERIC/fqZg3pt1YRvCKYL+bnEFbzoylxCaAWY3Uhz56oFRiQjsHuvJRGRtLRihbnBo60hJisxTGQQ39CU/FxiJ1WpxB5riMmf33ymQoxIxqAgIyLpUmwsjBwJbdva/p5ihBFKHQYxCYDJDKA2vxFGCfLlg1Gj4Nw5hRiRjCTF069FRFJD/FiYCRPg2jXb39ecVczjVfJwhcvk4VXmkbdLc+Y2gMKFNQ5GJKNSkBGRdCMkBHr2hIsXbX+PB9FM5C36Mh2ArTzFm/5LeXt6UbW8iGQCNgUZewb6Nm3aNMXFiEjmFRJijoWxR0mOsIx2VOF3ADZUfpss48fwW/2san0RySRsCjLNmzdPcGyxWLh7ZwOL5b9t7tPbFgUikv7Fxpqr9NqjDcv5nO54cZVLbvk5OHg+z455MXUKFJF0y6bBvnFxcdbXunXrqFixImvXriUyMpLIyEi+//57KleuzA8//JDa9YpIBhQaavtWA9m4yQxeZznt8OIqJ4vUxvvYXmoqxIhkSnaPkRkwYAAzZ86kdu3a1nMNGzYkR44c9OzZkwMHDji0QBHJ+CIibLuvDIdYTlsqsJ84LPzdYjCPLB8NWTTcTySzsnv69dGjR8mTJ0+i897e3hw/ftwBJYlIZlOgwIPv6cgidlOFCuznPL78NvQHHgkZqxAjksnZHWSqVavGm2++yblz56znzp07x9tvv82TTz7p0OJEJGOLjYXRo++/AWR2bjCH7iziZXJxndAs9dk1Zy9Pf/B82hUqIumW3f+UmTt3Li1atKBIkSIEBgYCcOrUKUqXLs3q1asdXZ+IZFC2TLUux18spy3l+ZM4LITWG07tH4fh7qEpSSJisjvIlCpViv379/PTTz9x8OBBAMqVK0eDBg0SzF4SEUmOLVOtgwjmU3qTkxtczVWQHCGLqPvcM2lToIi4DItx9zxqO926dQtPT890HWCioqLw9vbmypUreHl5ObsckUwvNhaKFUt+llIOrvMpvenCfAAuVm5A/u+/BD+/tCtSRJzO1p/fdo+RiYuL4/3336dw4cLkypWLsLAwAIYNG8YXX3yR8opFJEOLjYVNm6Bdu+RDTHn+xy6q0oX5xOLGe7zPT2/+oBAjIsmyO8h88MEHBAcH89FHH+Hh4WE9X758eT7//HOHFiciGUNIiNkKU78+rFyZ1B0G3ZnDDp6kHAcJpxD12cgY3qNgYY2HEZHk2R1kFixYwOzZs+nUqRPud60BXqFCBeuYGRGReCtWmONhkmuFyU0Ui+nIHHqSnVt8T2MqspdQniYw0NzsUUQkOXYHmfDwcEqVKpXofFxcHLdv33ZIUSKSMXz1FbRvn/z1CuxlF1XpwFLu4M7/8SFN+JZ/8MVigSlTtGO1iNyf3UHm0UcfJTQ0NNH5FStWUKlSJYcUJSKub8UKaNvWHBuTmMHrzGAbT1GGw5wkkKf5lY/4PwzcyJ/ffL92rxaRB7F7+vXw4cMJCgoiPDycuLg4QkJCOHToEAsWLODbb79NjRpFxMV89RV06JD0NS+u8DndacMKAL7mJboQzGXyAdClC3z+uVpiRMQ2drfINGvWjG+++Yb169eTM2dOhg8fzoEDB/jmm2947rnn7HpWsWLFsFgsiV69e/cGzOndvXv3Jn/+/OTKlYtWrVolWFFYRNKfkJDkW2KqsIs9VKINK4ghKwOZRDPWWENMQIBCjIjYx64WmTt37jB27Fi6du3KTz/99NAfvnPnTmLv+tvujz/+4LnnnqPNv+uVDxw4kO+++46vvvoKb29v+vTpQ8uWLdm8efNDf7aIOF5sLPTvn9QVg75MYyJv4cFtwihGO5axE3Nbk/ilqD75RCFGROxj94J4uXLl4o8//qBYsWIOL2bAgAF8++23HD58mKioKHx9fVm8eDGtW7cG4ODBg5QrV46tW7fy1FNP2fRMLYgnknY2bTKnWN8tD5eZS1dasBqAEFrQlblcIY/1nsBAc2CvxsSISLxUWxDv2Wef5Zdffnmo4pISExPDl19+SdeuXbFYLOzevZvbt2/ToEED6z2PPPIIRYoUYevWrQ7/fBF5eOHhCY+rs409VKIFq4nGgz5MoxUruUIe3NzM1puNGyEsTCFGRFLG7sG+jRs3ZvDgwfzvf/+jSpUq5MyZM8H1pk2bpqiQ1atXExkZSZcuXQA4e/YsHh4e5MmTJ8F9fn5+nD17NtnnREdHEx0dbT2OiopKUT0iYp+QEBgwwPy9hTgGMpkPGUxW7nCUErRlOb9TxXr/0qX33/VaRMQWdgeZN954A4BJkyYlumaxWBKMebHHF198QePGjSlUqFCK3h9v3LhxjBo16qGeISL2CQmB1q3BMCAfFwmmCy9hzmJcRlt6MpsovAFzDMzSpeb9IiIPK0V7LSX3SmmIOXHiBOvXr6d79+7WcwULFiQmJobIyMgE9547d46CBQsm+6whQ4Zw5coV6+vUqVMpqklEbBM/wNcwoBa/sZeKvMS33MKT15lBe5ZaQwzAkiUKMSLiOHYHmbvdunXLIUXMmzePAgUK8OKLL1rPValShaxZs7JhwwbruUOHDnHy5Elq1KiR7LM8PT3x8vJK8BKR1BMaCuGn4xjMODZRj0BOc4gyVGc7s3gdMKck+fqa+yypO0lEHMnuIBMbG5tg9+tjx44BKd/9Oi4ujnnz5hEUFESWLP/1dHl7e9OtWzfefPNNNm7cyO7du3n11VepUaOGzTOWRCT1XTp4nrU0ZhzvkoVYvqQTVdnFfiokuG/yZA3oFRHHszvIjBkzxqG7X69fv56TJ0/StWvXRNcmT55MkyZNaNWqFU8//TQFCxYkJCTE7s8QkVTyyy+8+F5FGrKOG2SnK1/wCgu5Ru5EtxYu7IT6RCTDs3sdmVKlSjFr1iyeffZZcufOzb59+yhRogQHDx6kRo0aXL58ObVqTRGtIyPiWLGxELoplrwzxvB4yCjcjDj+ohxtWc6flE90v8VirtgbFqbF7kTEdqm2jox2vxbJvEJCoFrgWWIbPE+FlSNwM+KYRxeqsTPJEBNPu1iLSGrR7tciYpMVK+CzVutZG1GBZ/mZ6+TgFRbQlXncIGeS73F3h+XLNTZGRFKPdr8WkQdasfQOhzqOYh1jcMNgP4/TjmUcpNx93xcbCz4+aVSkiGRKTt39WkTSt9hYmPxWOL4dnmWo8QFuGMymB9XZ/sAQEy8iIpWLFJFMze4WGYA6deo4ZPdrEUmfYmNhzBj434Qf+OzaK/jyD1fJRU9ms5QOdj3L3z+VihQR4SEXxBORjCM21ty9euBAKJD3NllHDOGra43x5R/2UJHK/G5XiLFYzF2t69RJvZpFRGxqkcmbNy8Wi8WmB166dOmhChKRtBcSYm4zcPo0BHCKNXSgNpsB+JQ3GMTHRJPN5ufF/3Wh2UoiktpsCjJTpkyx/v7ixYt88MEHNGzY0LpVwNatW/nxxx8ZNmxYqhQpIqnn7g0fX+Rb5hNEfi5xBS+68QUrsX9jpIAAM8RotpKIpDa7F8Rr1aoV9evXp0+fPgnOT58+nfXr17N69WpH1vfQtCCeSPJiYszQEXkhhnEMYRDmrvY7qUo7lhFGCZuf5esLnTpBs2Zmd5JaYkTkYdj689vuIJMrVy727t2baFG8I0eOULFiRa5du5ayilOJgoxI0kJC4LXXIOc/x1lKe55iOwBT6M//MZ4YPG16Tr9+0KKFwouIOFaqreybP39+1qxZk+j8mjVryJ8/v72PExEnWLECWrWCWv+sZg+VeIrtXCYPzVjNQKbYFGICAszdrD/5BOrVU4gREeewe/r1qFGj6N69O5s2baJ69eoAbN++nR9++IE5c+Y4vEARcYzYWAgNhVWrYPa0aKbwDv2ZCsA2qtOOZZykqE3PGjUKhg5VeBER57M7yHTp0oVy5coxdepU607U5cqV47fffrMGGxFJH+LDy5o1sGgRXLgAJThKKO2oym4AJvAW7zKWO2R94PPy54fZszWIV0TSD7uCzO3bt3nttdcYNmwYixYtSq2aRMQB7p5SHa81X/E53fEmiovkI4j5fEeTBz4rXz7zWWqFEZH0xq4xMlmzZmXlypWpVYuIOEj8lOr4EOPJLT7lDb6iLd5EEUptKrDvviHG2xsGDICNG+H8eRg+XCFGRNIfuwf7Nm/ePN1NsRaR/8TGmq0n8fMRS/M3W6nBG8wAYCxDqM9GwglI9hm+vmZ4mTxZA3lFJH2ze4xM6dKlGT16NJs3b6ZKlSrkzJkzwfV+/fo5rDgRsd+YMf+1xHRgMbN4jdxc4zy+vMJC1tHwgc+YORM8PFK5UBERB7B7HZnixYsn/zCLhWPHjj10UY6kdWQkMwkJMadVZ+cGn9CfHnwOwCbq0pHFRFDovu93d4elS81uKRERZ7L157fdLTJhYWEPVZiIpI74LqVHOMBy2vI4fxCHhfcZxmiGE8eD+4eWLFGIERHXkuLdr//55x/++ecfR9YiIg8hNBSeOT2fXVTlcf4ggoI8x0+MZNQDQ0xgoLm4XZs2aVSsiIiD2BVkIiMj6d27Nz4+Pvj5+eHn54ePjw99+vQhMjIylUoUkQe6do2iw4OYTxdycoOfaEBF9vIzzyb7Fl/f/2YlhYVpbRgRcU02dy1dunSJGjVqEB4eTqdOnShXrhwAf/31F8HBwWzYsIEtW7aQN2/eVCtWRJLwv/9B27YUP3iQWNwYyUjG8m6yrTCtW0Pv3tobSUQyBpsH+w4YMIANGzawfv16/Pz8Elw7e/Yszz//PM8++yyTJ09OlUJTSoN9JcMyDPj8c3PXxlu3MAoVou3txaz8py7J/b86IACOH1eAEZH0z+GbRq5evZqJEycmCjEABQsW5KOPPmLVqlUpq1ZE7HP1KnTqBD17wq1b0LAhcbv38nifpEOMxWK+PvlEIUZEMhabu5YiIiJ47LHHkr1evnx5zp4965CiROQ+9u6Ftm3h8GEMd3fCun7A9Bzv8OUTbly4kPRbAgJgyhSNgxGRjMfmIOPj48Px48cJCEh6NdCwsDDy5cvnsMJE5B6GYa5UN3AgREdzI38AndyWsnpOrfu+TTtVi0hGZnPXUsOGDRk6dCgxMTGJrkVHRzNs2DAaNWrk0OJE5F9XrkC7dvDGGxAdzd9lmhB4cS+rL9w/xFgs5jAaEZGMyubBvqdPn6Zq1ap4enrSu3dvHnnkEQzD4MCBA3z22WdER0eza9cuAgMDU7tmu2iwr7i8XbvMEHPsGHHuWXg/53hGRg0ELDY/YuNGc88kERFX4fCVfQMCAti6dStvvPEGQ4YMIT7/WCwWnnvuOaZPn57uQoyISzMMmDYN3noLbt/mum9Rnr2wjO1R1e1+VEREKtQnIpIO2LVFQfHixVm7di2XL1/m8OHDAJQqVUpjY0Qc7fJl6NoV/t1p3mjeguo7vuBPUrZOk7+/A2sTEUlH7N5rCSBv3rw8+eSTjq5FRAC2bze7kk6cMLegnjiRX8r34c9nbO9KimexmDOW6tRJhTpFRNKBFO+1JCIOZhjw8cdQu7YZYkqUgC1boG9fIs6mLMSAOe1aM5ZEJKNSkBFJDy5ehKZNzfEwd+6Y68T8/jtUqQKkrGsoIABWrNDaMSKSsaWoa0lEHGjzZmjfHk6fBk9Pswnltdf+a1LB7BoKCIDwcJLdfgDAxwdefhmaNdNeSiKSOahFRsRZ4uLgww+hbl0zxJQuDdu2weuvJwgxYAaSTz4xf29Jppdp1Cg4exYmTzanWivEiEhmoCAj4gznz8MLL8CQIRAbCx07wu7dULFism9p2dLsKipcOOH5wEBYuRKGD1d4EZHMR11LImntl1+gQwdzcZds2cy1Yrp1S76p5S4tW5rdRqGh5tv9/dWFJCKZm4KMSFqJjYWxY2HkSLNb6ZFH4KuvoHz5ZG+PDywFCpjnzp9XeBERuZuCjEhaOHvWHIW7YYN5HBQEn34KOXMmuC0+vKxZA4sWcd/drD/5RDOSREQUZERS24YN0KkTnDsHOXIQN/0zfi0eRMTXZutKzZrmcjEPCi93Cw+H1q01vVpEREFGJLXExppTiT74wJwzXb48P/VYTtfh5Th9+r/b3N3NW+1hGOaQmgEDzDEz6mYSkcxKQUYkNZw5Y85E+uUX87h7d9Y88wktOuVItA6MvSEmnmHAqVNmV5R2thaRzEpBRsTRfvgBXnkF/vkHcuWCWbOIad2RHgH3X8wupbSztYhkZlpHRsRR7twx14Vp3NgMMRUqwO7dhGTrSOHCto19SQntbC0imZmCjIgjnDpl9u98+KF53KsXbNtGyB9laN3azDWOZrGYi+FpZ2sRyczUtSTysL791pxOfekSeHnB559DmzbExkL//qnTnaSdrUVETGqREUmpmBhzt+qXXjJDTJUq5o7V/4aYadNIMDvJkbSztYiISS0yIilx/Li5Y/X27eZx//4wfjx4ehISYh46IsT4+ppL0DRpYh5rZV8RkYQUZETstXo1vPoqREZCnjwwbx6xLzW3rsg7ZYp9j3NzM3csiBcfXpo1U2AREXkQBRkRW0VHwzvvwNSpABhPVmf7wKUs+6UYi3qmbFaSry+cOGE27GgTSBER+ynIiNji2DFo2xZ27wbg76aDaLR7LGEdPFL0uPjBujNnQvbsWtBORCSlNNhX5EFWrIBKlcwQky8fW4Z8wyPfTCQsPGUhBjRYV0TEURRkRJJz6xa88Qa0aQNRUVCrFrG799JuYZOHmlI9eTKEhSnEiIg4goKMSFIOH4YaNWDGDPN4yBDYtInQ44Epno0Uv4Bd374aAyMi4igaIyNyryVLoGdPuHYNfHzgyy+hYUMg5fsaaQE7EZHUoRYZkXg3b0KPHuau1deuQd26sG8fNGxIbCxs2gR//ZWyR2tMjIhI6lCLjAjAgQPQrh38739m88mwYeYrS5aHWuBuwACtByMikpoUZETmzzcH9d64AX5+sGgRPPssACEh0Lq1/fslBQaa3UhqgRERSV1O71oKDw/n5ZdfJn/+/GTPnp3HH3+cXbt2Wa936dIFi8WS4NWoUSMnViwZxvXr0KWL+bpxwwwve/daQ4y9mz76+potMBs3alaSiEhacWqLzOXLl6lVqxb169dn7dq1+Pr6cvjwYfLmzZvgvkaNGjFv3jzrsaenZ1qXKhnNH3+YC9wdOGDuETByJLz7boL+n9BQ27qT3nvPzD7qPhIRSXtODTLjx48nMDAwQUgpXrx4ovs8PT0pWLBgWpYmGZVhwNy55hzomzehUCFYvNgc2HsPW2coPfqoVuYVEXEWp3Ytff3111StWpU2bdpQoEABKlWqxJw5cxLdt2nTJgoUKEDZsmXp1asXFy9eTPaZ0dHRREVFJXiJAHD1KrzyCnTvboaYhg3NrqQkQkxsLJw7Z9tj/f0dW6aIiNjOYhgPs0bpw8mWLRsAb775Jm3atGHnzp3079+fmTNnEhQUBMDSpUvJkSMHxYsX5+jRo7z77rvkypWLrVu34p5EO/7IkSMZNWpUovNXrlzBy8srdf9Akn7t22d2Jf39t9n/M2YMvP222a2EGVxCQ81WmMOHYc6cB3crWSzmtOqwMHUpiYg4WlRUFN7e3g/8+e3UIOPh4UHVqlXZsmWL9Vy/fv3YuXMnW7duTfI9x44do2TJkqxfv55n/x2Uebfo6Giio6Otx1FRUQQGBirIZFaGAbNmmaNwo6PN5LF0KdSqZb0lJdOr4xe409owIiKpw9Yg49SuJX9/fx599NEE58qVK8fJkyeTfU+JEiXw8fHhyJEjSV739PTEy8srwUsyqagoaN8eevUyQ0yTJmZX0j0hpnVr+9eI0QJ3IiLpg1ODTK1atTh06FCCc3///TdFixZN9j2nT5/m4sWL+GtggtzP7t1QuTIsXw5ZssDEifD115A/v/UWe6dXx9OmjyIi6YdTg8zAgQPZtm0bY8eO5ciRIyxevJjZs2fTu3dvAK5du8bbb7/Ntm3bOH78OBs2bKBZs2aUKlWKhv/ufSOSgGHAtGlQsyYcPQpFi8Jvv8GgQf/1B/3L1unV9/Lz05gYEZH0wqlBplq1aqxatYolS5ZQvnx53n//faZMmUKnTp0AcHd3Z//+/TRt2pQyZcrQrVs3qlSpQmhoqNaSkcQuX4ZWraBfP4iJgebNYc8eqF49wW3x+yatXJmyj1FjoIhI+uHUwb5pwdbBQuLiduww90o6fhyyZjW7kvr2TdQK8zD7JmmWkohI2rH157f2WhLXZhjmoJX/+z+4cwdKlIBly6Bq1US3pnTfJPgvD02ZohAjIpKeKMiI67p0ydwn6ZtvzOM2bcwFYLy9rbfErw8THg4DB6YsxIDZEqNNIEVE0h8FGXFNW7aYU6tPnQJPT7NV5vXXE3QlPUw3UkAA9OgBpUubY2K0j5KISPqkICOuJS4OJkyAoUPN5pbSpc0p1hUrJrgtpd1IffqY44UVXEREXIOCjLiOCxcgKAjWrjWPO3QwV+3NnTvBbSldHwbMEKMNIEVEXIeCjLiGX381g8uZM5Atm7lWTLduiWYlQcrWh4mfkVSnjoPqFRGRNOHUdWREHig2Fj74AOrXN0PMI4+YU627d08yxIC58aM9NCNJRMR1qUVG0q9z5+Dll2H9evM4KAg+/RRy5kx06927V587Z9/HaEaSiIjrUpCR9Onnn6FjRzOV5MgBn31mBpkkJDU7yd3dDDfJ8fU1JzoVLqyBvSIirkxBRtKX2FgYPRref98crfvYY+aspH93Sb+75cXfH/75B9q2TTywN7kQE9+NNHOmWmBERDICBRlJP86cgU6dzI2QwBwH88knZosMybe83G920r0tM+pGEhHJWBRkxCnubVmpc+NH3Lu8AhcuYOTKxYH+s9j3WEf8d5hdP2vWJL0uzP26j+KvT55s7lithe1ERDIeBRlxuEQh5Z7wcHfLijt3GM1w6jEOgMhiFWh6czmhY8pY7y9cGG7dSvn2An5+5sxtERHJeBRkxGYPCiiQdPdPQIDZQ9SyZcIVdwtzmiV0oA6/ATCD1xl4fDLRZEvwzPDwh6vb3//h3i8iIumXgozY5EEBJf6epLp/wsPN88uWwZtvmtdf4DvmE4QPF4kiN935nK9o69CatcidiEjGpwXx5IHiA8q9q+XGB5SQkPtvCxB/rndvOHv6Nh/xNt/RBB8uspvKVGJPqoQY0CJ3IiIZnVpk5L4eFFAsFhgwALy9778tgGFA9gsnWEN7arANgKn05W0mEIPnQ9ep2UkiIpmTgozc14P2LTIMOHXqvxnTyWnKGoLpQl4iuUweujKX1bR46PriW16WLDEXubvf+B0REcl4FGTkvuzdt+heWYnhI95hAJ8AsDvLk7S+s4zjFLPrORYL5MsH2bMnHqejlhcRkcxLY2Tkvmyd8VOvnhkq7t7HsTjH2Ewta4iZlXsQxxeEcsJSLNF+j3cfJ3dt9mw4fhw2boTFi81fw8IUYkREMjMFGbmvOnUSB5S7WSwQGGgGmU8++e9cS1ayh0pUYxcXyUdTvsY3eCKtOniwYoW5NszdAgJg5UrzldS1FSvMwOLubn5Whw7mr+o+EhHJ3CyGkdJlxlxDVFQU3t7eXLlyBS8vL2eX45LiZy1BwkG/8eEmPmQArF56i8jub9Hl+qcA/EYtBvkv4f+mByZoObnfmjS2rFcjIiIZm60/vxVkxCZJrSMTGHjP+JTDh6FdO9izB4C/mg7mQt/R1K6fVUFERETsoiDzLwUZx7lvS8nSpdCjB1y7Bj4+sHAhNGrk1HpFRMR12frzW7OWMqmUdN/Ej09J4OZNcyGZ2bPN46efNkfi3jvQRUREJBVosG8mFBICxYpB/frQsaP5a7Fi5nm7HDwI1aubIcZigWHDYMMGhRgREUkzCjKZjC3bDdhk4UKoWhX+9z9ze+l162D0aMiiRj4REUk7CjKZiC37IQ0YkHCp/0SuX4dXX4XOnc3fP/MM7N0LDRqkQsUiIiL3pyCTidi63UBoaDI3/PknPPkkBAeDm5vZArNuHRQsmBrlioiIPJD6AVzMw6yxYut2A4nuMwyYNw/69DEH9/r7mwN6E438FRERSVsKMi4kqbVcAgLMFXVtWabf1u0GEtx37Rq8/josWmQeP/+8OT6mQAGb6xYREUkt6lpyEY4YpGvrdgN16vx7Yt8+qFLFDDHu7jBuHKxdqxAjIiLphoKMC3DIIF3MLHL3fkh3iz+eMgXc3QyYNcucWv3332b62bQJBg82x8aIiIikE/qp5AIeepDuXVq2JNlNG1esgJYNoswdGV9/HaKj4cUXzS0Hatd+uD+EiIhIKtAYGReQ4kG6yWjZEpo1S2LQ8L7foXJbOHrUXA9m3Dh48021woiISLqlIOMCUjRI9wESbDdgGPDppzBoEMTEQNGi5t5JTz1lb6kiIiJpSv/UdgF2D9K1R2SkOVq4b18zxDRvbnYlKcSIiIgLUJBxATYP0rVxPRmrHTugUiVzylPWrOaHhIRA3rwPW7KIiEiaUJBxEQ8cpGvDOjJWhgGTJ5sDeI8fhxIlYMsW6Ncv+WYfERGRdEhjZFxIsoN07WmJuXQJunSBb74xj1u3hs8/B2/v1ChZREQkVSnIuJgEg3TttWULtG9vztX29DRbZV5/Xa0wIiListS1lBnExcFHH8HTT5shpnRp2LYNevVSiBEREZemFpkUeJiNG9PchQsQFGRuLQDmYnezZkHu3M6tS0RExAEUZOz0sBs3pqlffzWDy5kzkC0bTJ0K3burFUZERDIMdS3ZwREbN6aJuDgYMwbq1zdDTNmysH079OihECMiIhmKgoyNHLVxY6o7dw4aNYL33jMDzSuvwK5d8MQTTi5MRETE8RRkbOTIjRtTzc8/Q8WK8NNPkD07zJsHCxZArlxOLEpERCT1KMjYyNEbNzpUbCyMHAkNGsDZs/DYY2YrTJcuTihGREQk7Wiwr41SY+NGh4iIgI4dYdMm87hrV5g2DXLkSONCRERE0p5aZGyUqhs3ptS6dVChghlicuaEhQvhiy8UYkREJNNQkLFRqm3cmBJ37sDQoeag3gsXzIG8u3fDyy+nwYeLiIikHwoydnDoxo0pdfo0PPMMjB1rjjB+7TVzld6yZdPgw0VERNIXjZGxk0M2bkyp77+Hzp3h4kVzZd45c6BduzT4YBERkfRJQSYFHmrjxpS4fdvsSpowwTyuXBmWLYNSpdKwCBERkfRHQSa9O3HC3LF62zbzuE8fmDjR3L1aREQkk1OQSc/WrIFXX4XLl8Hb25yR1KqVs6sSERFJNzTYNz2KiYGBA6F5czPEVKsGe/YoxIiIiNxDQSa9CQuD2rXNudxgBprffoPixZ1aloiISHrk9CATHh7Oyy+/TP78+cmePTuPP/44u3btsl43DIPhw4fj7+9P9uzZadCgAYcPH3Zixalo5UqoVAl27oS8ec2upUmTwMPD2ZWJiIikS04NMpcvX6ZWrVpkzZqVtWvX8tdff/Hxxx+TN29e6z0fffQRU6dOZebMmWzfvp2cOXPSsGFDbt265cTKHezWLXMQb+vWcOUK1KgBe/dC06bOrkxERCRdsxiGYTjrwwcPHszmzZsJTWbLaMMwKFSoEIMGDeKtt94C4MqVK/j5+REcHEz79u0f+BlRUVF4e3tz5coVvLy8HFq/Qxw5Am3bmmNgAN55Bz74ALJmdW5dIiIiTmTrz2+ntsh8/fXXVK1alTZt2lCgQAEqVarEnDlzrNfDwsI4e/YsDRo0sJ7z9vamevXqbN26NclnRkdHExUVleCVbi1daq4Js2cP+PiYC96NH68QIyIiYiOnBpljx44xY8YMSpcuzY8//kivXr3o168f8+fPB+Ds2bMA+Pn5JXifn5+f9dq9xo0bh7e3t/UVGBiYun+IlLh509xaoEMHuHrVXBp4715o3NjZlYmIiLgUpwaZuLg4KleuzNixY6lUqRI9e/akR48ezJw5M8XPHDJkCFeuXLG+Tp065cCKHeDQIXjqKZg929xt8r334OefE2/gJCIiIg/k1CDj7+/Po48+muBcuXLlOHnyJAAFCxYE4Ny5cwnuOXfunPXavTw9PfHy8krwSje+/BKqVIH9+6FAAfjxR3j/fciidQlFRERSwqlBplatWhw6dCjBub///puiRYsCULx4cQoWLMiGDRus16Oioti+fTs1atRI01ofyo0b0LUrvPIKXL8O9eubXUnPPefsykRERFyaU5sCBg4cSM2aNRk7dixt27Zlx44dzJ49m9mzZwNgsVgYMGAAH3zwAaVLl6Z48eIMGzaMQoUK0bx5c2eWbrs//zRnJf31F7i5wYgR5gaQabJdtoiISMbm1CBTrVo1Vq1axZAhQxg9ejTFixdnypQpdOrUyXrPO++8w/Xr1+nZsyeRkZHUrl2bH374gWzZsjmxchsYBgQHQ+/e5uDeggVhyZI03jZbREQkY3PqOjJpwSnryFy7Bm+8AQsXmsfPP2/+vkCBtPl8ERERF+cS68hkSPv3Q9WqZnBxd4exY2HtWoUYERGRVKDpMo5iGDBnDvTvb245ULiwueBd7drOrkxERCTDUpBxhKgoc4G7pUvN4xdegPnzzdV6RUREJNWoa+lh7dljrg2zdKm5HsyECfDNNwoxIiIiaUAtMillGPDZZ/DmmxATA0WKmGHGlda3ERERcXEKMilhGObidosWmcfNmsHcuZAvn3PrEhERyWTUtZQSFou5X1LWrDBlCqxapRAjIiLiBGqRSanevc31YcqUcXYlIiIimZZaZFLKYlGIERERcTIFGREREXFZCjIiIiLishRkRERExGUpyIiIiIjLUpARERERl6UgIyIiIi5LQUZERERcloKMiIiIuCwFGREREXFZCjIiIiLishRkRERExGUpyIiIiIjLUpARERERl5XF2QWkNsMwAIiKinJyJSIiImKr+J/b8T/Hk5Phg8zVq1cBCAwMdHIlIiIiYq+rV6/i7e2d7HWL8aCo4+Li4uI4c+YMuXPnxmKxOLscp4iKiiIwMJBTp07h5eXl7HIyNH3XaUvfd9rRd5229H2bLTFXr16lUKFCuLklPxImw7fIuLm5ERAQ4Owy0gUvL69M+3+ItKbvOm3p+047+q7TVmb/vu/XEhNPg31FRETEZSnIiIiIiMtSkMkEPD09GTFiBJ6ens4uJcPTd5229H2nHX3XaUvft+0y/GBfERERybjUIiMiIiIuS0FGREREXJaCjIiIiLgsBRkRERFxWQoymcjx48fp1q0bxYsXJ3v27JQsWZIRI0YQExPj7NIyrDFjxlCzZk1y5MhBnjx5nF1OhvLpp59SrFgxsmXLRvXq1dmxY4ezS8qwfv31V1566SUKFSqExWJh9erVzi4pwxo3bhzVqlUjd+7cFChQgObNm3Po0CFnl5WuKchkIgcPHiQuLo5Zs2bx559/MnnyZGbOnMm7777r7NIyrJiYGNq0aUOvXr2cXUqGsmzZMt58801GjBjB77//ToUKFWjYsCHnz593dmkZ0vXr16lQoQKffvqps0vJ8H755Rd69+7Ntm3b+Omnn7h9+zbPP/88169fd3Zp6ZamX2dyEyZMYMaMGRw7dszZpWRowcHBDBgwgMjISGeXkiFUr16datWqMX36dMDcUy0wMJC+ffsyePBgJ1eXsVksFlatWkXz5s2dXUqmcOHCBQoUKMAvv/zC008/7exy0iW1yGRyV65cIV++fM4uQ8RmMTEx7N69mwYNGljPubm50aBBA7Zu3erEykQc78qVKwD6e/o+FGQysSNHjjBt2jRee+01Z5ciYrN//vmH2NhY/Pz8Epz38/Pj7NmzTqpKxPHi4uIYMGAAtWrVonz58s4uJ91SkMkABg8ejMViue/r4MGDCd4THh5Oo0aNaNOmDT169HBS5a4pJd+3iIi9evfuzR9//MHSpUudXUq6lsXZBcjDGzRoEF26dLnvPSVKlLD+/syZM9SvX5+aNWsye/bsVK4u47H3+xbH8vHxwd3dnXPnziU4f+7cOQoWLOikqkQcq0+fPnz77bf8+uuvBAQEOLucdE1BJgPw9fXF19fXpnvDw8OpX78+VapUYd68ebi5qVHOXvZ83+J4Hh4eVKlShQ0bNlgHnMbFxbFhwwb69Onj3OJEHpJhGPTt25dVq1axadMmihcv7uyS0j0FmUwkPDycevXqUbRoUSZOnMiFCxes1/Qv2dRx8uRJLl26xMmTJ4mNjWXv3r0AlCpVily5cjm3OBf25ptvEhQURNWqVXnyySeZMmUK169f59VXX3V2aRnStWvXOHLkiPU4LCyMvXv3ki9fPooUKeLEyjKe3r17s3jxYtasWUPu3Lmt4768vb3Jnj27k6tLpwzJNObNm2cASb4kdQQFBSX5fW/cuNHZpbm8adOmGUWKFDE8PDyMJ5980ti2bZuzS8qwNm7cmOT/joOCgpxdWoaT3N/R8+bNc3Zp6ZbWkRERERGXpQESIiIi4rIUZERERMRlKciIiIiIy1KQEREREZelICMiIiIuS0FGREREXJaCjIiIiLgsBRkRF7dp0yYsFguRkZHOLsUuFouF1atXO+x5xYoVY8qUKQ57Xlo7fvw4FovFuvqzq/53FUlrCjIi6diDdtkeOXKks0t8oJEjR1KxYsVE5yMiImjcuHHaF5QOdOnSxbpPVLzAwEAiIiIoX768c4oScVHaa0kkHYuIiLD+ftmyZQwfPpxDhw5Zz+XKlYtdu3Y5ozRiYmLw8PBI8fu1v1dC7u7u+k5EUkAtMiLpWMGCBa0vb29vLBZLgnN3bzy5e/duqlatSo4cOahZs2aCwAOwZs0aKleuTLZs2ShRogSjRo3izp071usnT56kWbNm5MqVCy8vL9q2bcu5c+es1+NbVj7//HOKFy9OtmzZAIiMjKR79+74+vri5eXFM888w759+wAIDg5m1KhR7Nu3z9qKFBwcDCTuWjp9+jQdOnQgX7585MyZk6pVq7J9+3YAjh49SrNmzfDz8yNXrlxUq1aN9evX2/VdxsbG8uabb5InTx7y58/PO++8Q1BQUIKWkaS6pypWrJig5WvSpEk8/vjj5MyZk8DAQN544w2uXbtmvR4cHEyePHn48ccfKVeuHLly5aJRo0bWUDpy5Ejmz5/PmjVrrN/Jpk2bEnUtJeW3336jTp06ZM+encDAQPr168f169et1z/77DNKly5NtmzZ8PPzo3Xr1nZ9RyKuSEFGJIMYOnQoH3/8Mbt27SJLlix07drVei00NJTOnTvTv39//vrrL2bNmkVwcDBjxowBIC4ujmbNmnHp0iV++eUXfvrpJ44dO0a7du0SfMaRI0dYuXIlISEh1h+4bdq04fz586xdu5bdu3dTuXJlnn32WS5dukS7du0YNGgQjz32GBEREURERCR6Jpi7K9etW5fw8HC+/vpr9u3bxzvvvENcXJz1+gsvvMCGDRvYs2cPjRo14qWXXuLkyZM2fz8ff/wxwcHBzJ07l99++41Lly6xatUqe79m3NzcmDp1Kn/++Sfz58/n559/5p133klwz40bN5g4cSILFy7k119/5eTJk7z11lsAvPXWW7Rt29YabiIiIqhZs+YDP/fo0aM0atSIVq1asX//fpYtW8Zvv/1Gnz59ANi1axf9+vVj9OjRHDp0iB9++IGnn37a7j+fiMtx9q6VImKbefPmGd7e3onOx+9MvH79euu57777zgCMmzdvGoZhGM8++6wxduzYBO9buHCh4e/vbxiGYaxbt85wd3c3Tp48ab3+559/GoCxY8cOwzAMY8SIEUbWrFmN8+fPW+8JDQ01vLy8jFu3biV4dsmSJY1Zs2ZZ31ehQoVEdQPGqlWrDMMwjFmzZhm5c+c2Ll68aOO3YRiPPfaYMW3aNOtx0aJFjcmTJyd7v7+/v/HRRx9Zj2/fvm0EBAQYzZo1u+8zKlSoYIwYMSLZ53711VdG/vz5rcfxu8wfOXLEeu7TTz81/Pz8rMdBQUEJPtcwDCMsLMwAjD179hiG8d9/18uXLxuGYRjdunUzevbsmeA9oaGhhpubm3Hz5k1j5cqVhpeXlxEVFZVsrSIZkcbIiGQQTzzxhPX3/v7+AJw/f54iRYqwb98+Nm/ebG2BAbOr5datW9y4cYMDBw4QGBhIYGCg9fqjjz5Knjx5OHDgANWqVQOgaNGi+Pr6Wu/Zt28f165dI3/+/AlquXnzJkePHrW59r1791KpUiXy5cuX5PVr164xcuRIvvvuOyIiIrhz5w43b960uUXmypUrREREUL16deu5LFmyULVqVQzDsLlOgPXr1zNu3DgOHjxIVFQUd+7csX6POXLkACBHjhyULFnS+h5/f3/Onz9v1+fca9++fezfv59FixZZzxmGQVxcHGFhYTz33HMULVqUEiVK0KhRIxo1akSLFi2sNYlkVAoyIhlE1qxZrb+3WCwACbpmRo0aRcuWLRO9L36siy1y5syZ4PjatWv4+/uzadOmRPfmyZPH5udmz579vtffeustfvrpJyZOnEipUqXInj07rVu3JiYmxubPsIWbm1uiYHP79m3r748fP06TJk3o1asXY8aMIV++fPz2229069aNmJgYa2i4+78FmP897A1M97p27RqvvfYa/fr1S3StSJEieHh48Pvvv7Np0ybWrVvH8OHDGTlyJDt37rTrv4WIq1GQEckEKleuzKFDhyhVqlSS18uVK8epU6c4deqUtVXmr7/+IjIykkcfffS+zz179ixZsmShWLFiSd7j4eFBbGzsfet74okn+Pzzz7l06VKSrTKbN2+mS5cutGjRAjB/qB8/fvy+z7ybt7c3/v7+bN++3Tpu5M6dO9YxPfF8fX0TzBSLiooiLCzMerx7927i4uL4+OOPcXMzhxguX77c5jri2fKd3Kty5cr89ddfyf43BLOVqUGDBjRo0IARI0aQJ08efv755yQDrEhGocG+IpnA8OHDWbBgAaNGjeLPP//kwIEDLF26lPfeew+ABg0a8Pjjj9OpUyd+//13duzYQefOnalbty5Vq1ZN9rkNGjSgRo0aNG/enHXr1nH8+HG2bNnC0KFDrdPCixUrRlhYGHv37uWff/4hOjo60XM6dOhAwYIFad68OZs3b+bYsWOsXLmSrVu3AlC6dGnrAON9+/bRsWNHa2uTrfr378+HH37I6tWrOXjwIG+88UaixeaeeeYZFi5cSGhoKP/73/8ICgrC3d3der1UqVLcvn2badOmcezYMRYuXMjMmTPtqiP+O9m/fz+HDh3in3/+SdDqk5z/+7//Y8uWLfTp04e9e/dy+PBh1qxZYx3s++233zJ16lT27t3LiRMnWLBgAXFxcZQtW9bu+kRciYKMSCbQsGFDvv32W9atW0e1atV46qmnmDx5MkWLFgXMro81a9aQN29enn76aRo0aECJEiVYtmzZfZ9rsVj4/vvvefrpp3n11VcpU6YM7du358SJE/j5+QHQqlUrGjVqRP369fH19WXJkiWJnuPh4cG6desoUKAAL7zwAo8//jgffvihNURMmjSJvHnzUrNmTV566SUaNmyYoCXFFoMGDeKVV14hKCiIGjVqkDt3bmsLT7whQ4ZQt25dmjRpwosvvkjz5s0TjHWpUKECkyZNYvz48ZQvX55FixYxbtw4u+oA6NGjB2XLlqVq1ar4+vqyefPmB77niSee4JdffuHvv/+mTp06VKpUieHDh1OoUCHA7MoLCQnhmWeeoVy5csycOZMlS5bw2GOP2V2fiCuxGA/bcSsi4qK6dOlCZGSkQ7dKEJG0pRYZERERcVkKMiIiIuKy1LUkIiIiLkstMiIiIuKyFGRERETEZSnIiIiIiMtSkBERERGXpSAjIiIiLktBRkRERFyWgoyIiIi4LAUZERERcVkKMiIiIuKy/h/ycSSGOcT8owAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "f = stats.probplot(smpl, dist = \"norm\", plot = pylab)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Чтобы построить QQ-график, нужно:\n",
    "\n",
    "* отсортировать выборку по неубыванию.\n",
    "* каждая точка выборки ставится на QQ-график, значения отображаются по оси Y\n",
    "* значение по Х - мат.ожидание квантиля распределения N(0,1), посчитанного по выборке такого как в нашем эксперименте объема.\n",
    "\n",
    "Понятно, что в случае нормальной выборки этот график должен лечь на прямую или, по крайней мере, быть близок к ней. Если ваши точки лучше описываются не прямой, а какой-то другой функцией, то наверняка ваше выборка далека от нормальности. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Оцифровать отклонения выборки от QQ-графика можно с помощью критерия Шапиро-Уилка. Выглядит он так:\n",
    "\n",
    "## Критерий Шапиро-Уилка\n",
    "\n",
    "\n",
    "**Выборка**    $X_n = {X_1 ... X_n}$ \n",
    "\n",
    "**Нулевая гипотеза**    $H_0: X \\sim N(\\mu,\\sigma^2)$\n",
    "\n",
    "**Альтернативная гипотеза**     $H_1: X \\not\\sim N(\\mu,\\sigma^2)$ \n",
    "\n",
    "Превращаем нашу выборку в вариационный ряд и строим QQ-график.\n",
    "\n",
    "**Статистика** $W(X_n) = \\dfrac{(\\sum^{n}_{i=1}a_iX_{(i)})^2}{\\sum^{n}_{i=1}(X_i - \\overline{X})^2}$\n",
    "\n",
    "$a_i$ - некие табулированые величины, связанные с матожиданиями квантилей $N(0,1)$\n",
    "\n",
    "**Нулевое распределение** Табличное, в аналитическом виде не выражается.\n",
    "\n",
    "Несмотря на то, что аналитического выражения нет, и распределение и соответствующие p-value табулированы, поэтому применять такой критерий нам ничто не помешает."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Давайте проверим нормальность с помощью критерия Шапиро-Уилка."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "ShapiroResult(statistic=0.9829477667808533, pvalue=0.22395388782024384)"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stats.shapiro(smpl)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Нормальность снова не отвергается :)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## критерии различаются по мощности"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Критерий Jarque-Bera\n",
    "Тест Жарка-Бера проверяет, соответствуют ли выборочные данные по асимметрии и эксцессу нормальному распределению.\n",
    "Следуя документации, устойчивые результаты этот критерий показывает с выборками объема от 2000 элеменетов. Но если в выборке сильно нарушены асимметрия и эксцесс, то можно заметить разницу в выборках более малого объема.\n",
    "\n",
    "**Выборка**    $X_n = {X_1 ... X_n}$ \n",
    "\n",
    "**Нулевая гипотеза**    $H_0: X \\sim N(\\mu,\\sigma^2)$\n",
    "\n",
    "**Альтернативная гипотеза**     $H_1: X \\not\\sim N(\\mu,\\sigma^2)$ \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "SignificanceResult(statistic=2.613253660022396, pvalue=0.27073174207220996)"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stats.jarque_bera(smpl)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "smpl2 = np.array([pylab.gamma(70, 5) for _ in range(500)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "SignificanceResult(statistic=2.1964466376078624, pvalue=0.33346301517425514)"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stats.jarque_bera(smpl2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## А зачем нам вообще проверять нормальность?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Часто используются тесты, предполагающие нормальность данных, но t-тесты, например, работают и в случае некоторого отклонения от нормальности. В реальной жизни нормальных выборок не существует, потому что нормальное распределение - модель, она должна **в каком-то смысле** адекватно описывать реальную ситуацию. Кроме того, проверка нормальности - тоже статистический тест, а значит на маленьких выборках мы почти никогда не сможем отвергнуть нормальность, а на очень больших выборках – почти всегда будем ее отвергать (оцени критерий Жарки-Бера). Поэтому предлагаю следующий способ использования:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1) Смотрим на специфику наших данных. Если они категориальные или бинарные, лучше использовать специальные критерии для них.\n",
    "\n",
    "2) Если величина непрерывная, постройте QQ-график. Если на нем выборка не очень сильно отклоняется от прямой, можно смело применять устойчивые к небольшим отклонениям методы. Например критерий Стьюдента.\n",
    "\n",
    "3) Если ваш метод все-таки чувствителен к нормальности, оцифруйте отклонение с помощью критерия Шапиро-Уилка \n",
    "\n",
    "4) Если критерий Шапиро-Уилка отвергает нормальность, лучше применить другой метод (например непараметрический, который не делает никаких предположений о конкретном характере распределения)\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}